.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "api/gallery/signal_processing/spectral/plot_07_yule_autoregression.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code or to run this example in your browser via JupyterLite. .. rst-class:: sphx-glr-example-title .. _sphx_glr_api_gallery_signal_processing_spectral_plot_07_yule_autoregression.py: Yule's autoregressive model (1927) ================================== Schuster's periodogram looked for strict periodicities hidden in noise. George Udny Yule, studying Wolfer's sunspot numbers, proposed another picture: a pendulum "pelted with peas", a damped oscillator kept going by random shocks, :math:`x_n = -a_1 x_{n-1} - a_2 x_{n-2} + e_n`. Its cycle wanders in phase and amplitude, as the sunspot cycle does, and it is fitted from the record's correlations, not its Fourier coefficients. Gilbert Walker (1931) extended the equations to any order. The fitted model gives a smooth spectrum :math:`\sigma^2/|A(e^{j\omega})|^2` (:func:`~engineeringkit.signal_processing.yule_walker`, :class:`~engineeringkit.signal_processing.ARModel`). .. GENERATED FROM PYTHON SOURCE LINES 18-24 .. code-block:: Python import matplotlib.pyplot as plt import numpy as np from scipy import signal from engineeringkit.signal_processing import power_spectral_density, yule_walker .. GENERATED FROM PYTHON SOURCE LINES 25-30 A disturbed pendulum with an 11-year period ------------------------------------------- Poles at radius 0.9 and angle :math:`2\pi/11` per year give a damped oscillator; white shocks keep it ringing. 300 years of annual values. .. GENERATED FROM PYTHON SOURCE LINES 30-40 .. code-block:: Python radius, period = 0.9, 11.0 a_true = np.array([1.0, -2 * radius * np.cos(2 * np.pi / period), radius**2]) rng = np.random.default_rng(1927) x = signal.lfilter([1.0], a_true, rng.normal(size=600))[300:] + 50.0 years = 1700 + np.arange(x.size) fig, axes = plt.subplots(1, 2, figsize=(11, 4)) axes[0].plot(years, x) axes[0].set(xlabel="year", ylabel="activity", title="Yule's process: no fixed phase or amplitude") .. image-sg:: /api/gallery/signal_processing/spectral/images/sphx_glr_plot_07_yule_autoregression_001.png :alt: Yule's process: no fixed phase or amplitude :srcset: /api/gallery/signal_processing/spectral/images/sphx_glr_plot_07_yule_autoregression_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none [Text(0.5, 14.833333333333327, 'year'), Text(104.22222222222221, 0.5, 'activity'), Text(0.5, 1.0, "Yule's process: no fixed phase or amplitude")] .. GENERATED FROM PYTHON SOURCE LINES 41-43 Fit by the Yule-Walker equations -------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 43-51 .. code-block:: Python model = yule_walker(x, 2) print("fitted A(z):", np.round(model.coefficients, 3), " true:", np.round(a_true, 3)) print(model.summary()) poles = np.roots(model.coefficients) fitted_period = 2 * np.pi / np.abs(np.angle(poles[0])) print(f"pole radius {abs(poles[0]):.3f}, period {fitted_period:.1f} years") .. rst-class:: sphx-glr-script-out .. code-block:: none fitted A(z): [ 1. -1.516 0.812] true: [ 1. -1.514 0.81 ] {'order': 2, 'noise_variance': 1.0905233922949178, 'prediction_gain': 9.765797731115361, 'stable': True} pole radius 0.901, period 11.0 years .. GENERATED FROM PYTHON SOURCE LINES 52-57 AR spectrum against the periodogram ----------------------------------- The periodogram of 300 values is ragged at every frequency; the two-parameter model gives a smooth peak at the same place. .. GENERATED FROM PYTHON SOURCE LINES 57-73 .. code-block:: Python f, P = power_spectral_density(x, 1.0) f_ar, P_ar = model.psd(fs=1.0, n=1024) true_psd = 2 / np.abs(np.exp(-2j * np.pi * np.outer(f_ar, np.arange(3))) @ a_true) ** 2 axes[1].semilogy(f[1:], P[1:], color="0.7", lw=0.7, label="periodogram") axes[1].semilogy(f_ar, P_ar, label="AR(2) fit") axes[1].semilogy(f_ar, true_psd, "k--", lw=0.8, label="true spectrum") axes[1].set(xlabel="frequency [1/year]", ylabel="PSD [units²·year]", title="Spectrum") axes[1].legend() fig.tight_layout() peak_period = 1 / f_ar[np.argmax(P_ar)] print(f"AR spectrum peaks at a period of {peak_period:.1f} years") assert np.allclose(model.coefficients, a_true, atol=0.06) assert abs(fitted_period - period) < 1.0 and abs(peak_period - period) < 1.5 plt.show() .. rst-class:: sphx-glr-script-out .. code-block:: none AR spectrum peaks at a period of 11.2 years .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.200 seconds) .. _sphx_glr_download_api_gallery_signal_processing_spectral_plot_07_yule_autoregression.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: lite-badge .. image:: images/jupyterlite_badge_logo.svg :target: ../../../../lite/lab/index.html?path=api/gallery/signal_processing/spectral/plot_07_yule_autoregression.ipynb :alt: Launch JupyterLite :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_07_yule_autoregression.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_07_yule_autoregression.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_07_yule_autoregression.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_