.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "api/gallery/statistics/time_series/plot_01_yule_autoregression.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_api_gallery_statistics_time_series_plot_01_yule_autoregression.py>`
        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_statistics_time_series_plot_01_yule_autoregression.py:


Yule's autoregressive model of the sunspot cycle
================================================

Yule (1927) argued that Wolfer's sunspot numbers were not a clean sine
wave hidden in noise, but a pendulum kicked at random. He fitted

.. math::

   X_t = \phi_1 X_{t-1} + \phi_2 X_{t-2} + \varepsilon_t,
   \qquad \phi_1 \approx 1.34,\ \phi_2 \approx -0.65,

the first autoregressive model, and found its period from the roots of
:math:`z^2 - \phi_1 z - \phi_2`. Here a series is simulated from Yule's
coefficients (no sunspot data ship with mathematicskit), and the
Yule-Walker equations recover them from the series alone.

.. GENERATED FROM PYTHON SOURCE LINES 20-29

.. code-block:: Python

    import matplotlib.pyplot as plt
    import numpy as np

    from mathematicskit.statistics import ar_autocorrelation, autocorrelation, simulate_ar, yule_walker

    phi = np.array([1.34, -0.65])
    years = 300
    series = 50 + simulate_ar(phi, n=years, sigma=15.0, seed=1927)








.. GENERATED FROM PYTHON SOURCE LINES 30-35

The implied cycle
-------------------

The characteristic roots are complex, :math:`r e^{\pm i\theta}`: the
model is a damped oscillation with period :math:`2\pi/\theta`.

.. GENERATED FROM PYTHON SOURCE LINES 35-40

.. code-block:: Python


    roots = np.roots([1.0, -phi[0], -phi[1]])
    period = 2 * np.pi / abs(np.angle(roots[0]))
    print(f"roots {roots.round(3)}, modulus {abs(roots[0]):.3f}, period {period:.1f} years (sunspot cycle: about 11)")





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    roots [0.67+0.448j 0.67-0.448j], modulus 0.806, period 10.7 years (sunspot cycle: about 11)




.. GENERATED FROM PYTHON SOURCE LINES 41-43

The Yule-Walker fit
---------------------

.. GENERATED FROM PYTHON SOURCE LINES 43-51

.. code-block:: Python


    fit = yule_walker(series, 2)
    fit_period = 2 * np.pi / abs(np.angle(np.roots([1.0, -fit.coefficients[0], -fit.coefficients[1]])[0]))
    print(f"fitted phi = {fit.coefficients.round(3)} (true {phi}), noise sd {np.sqrt(fit.noise_variance):.1f} (true 15.0)")
    print(f"fitted period {fit_period:.1f} years")
    for order in (1, 2, 3, 4):
        print(f"AR({order}) innovation variance: {yule_walker(series, order).noise_variance:7.1f}")





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    fitted phi = [ 1.362 -0.666] (true [ 1.34 -0.65]), noise sd 15.6 (true 15.0)
    fitted period 10.8 years
    AR(1) innovation variance:   440.2
    AR(2) innovation variance:   244.8
    AR(3) innovation variance:   244.4
    AR(4) innovation variance:   244.2




.. GENERATED FROM PYTHON SOURCE LINES 52-54

Series and autocorrelation
----------------------------

.. GENERATED FROM PYTHON SOURCE LINES 54-71

.. code-block:: Python


    fig, (ax0, ax1) = plt.subplots(1, 2, figsize=(14, 4.2), gridspec_kw={"width_ratios": [1.6, 1]})
    ax0.plot(np.arange(years), series, color="tab:orange")
    ax0.set_xlabel("year")
    ax0.set_ylabel("simulated sunspot number")
    ax0.set_title("A pendulum kicked at random: AR(2) with Yule's coefficients")

    lags = np.arange(41)
    ax1.stem(lags, autocorrelation(series, 40), basefmt=" ", label="sample ACF")
    ax1.plot(lags, ar_autocorrelation(phi, 40), "k-", label="AR(2) theory")
    ax1.axhspan(-2 / np.sqrt(years), 2 / np.sqrt(years), color="0.85", label=r"$\pm 2/\sqrt{n}$")
    ax1.set_xlabel("lag (years)")
    ax1.set_title(f"Damped oscillation, period {period:.1f} years")
    ax1.legend(fontsize=8)
    fig.tight_layout()

    plt.show()



.. image-sg:: /api/gallery/statistics/time_series/images/sphx_glr_plot_01_yule_autoregression_001.png
   :alt: A pendulum kicked at random: AR(2) with Yule's coefficients, Damped oscillation, period 10.7 years
   :srcset: /api/gallery/statistics/time_series/images/sphx_glr_plot_01_yule_autoregression_001.png
   :class: sphx-glr-single-img






.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 0.077 seconds)


.. _sphx_glr_download_api_gallery_statistics_time_series_plot_01_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/statistics/time_series/plot_01_yule_autoregression.ipynb
        :alt: Launch JupyterLite
        :width: 150 px

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: plot_01_yule_autoregression.ipynb <plot_01_yule_autoregression.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: plot_01_yule_autoregression.py <plot_01_yule_autoregression.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: plot_01_yule_autoregression.zip <plot_01_yule_autoregression.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_