.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_data_analysis/estimate_stochastic_processes/plot_estimate_arma.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 .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_data_analysis_estimate_stochastic_processes_plot_estimate_arma.py: Estimate a scalar ARMA process ============================== .. GENERATED FROM PYTHON SOURCE LINES 6-54 The objective here is to estimate an ARMA model from a scalar stationary time series using the Whittle estimator and a centered normal white noise. The data can be a unique time series or several time series collected in a process sample. If the user specifies the order :math:`(p,q)`, OpenTURNS fits a model :math:`ARMA(p,q)` to the data by estimating the coefficients :math:`(a_1, \dots, a_p), (b_1, \dots, b_q)` and the variance :math:`\sigma` of the white noise. If the User specifies a range of orders :math:`(p,q) \in Ind_p \times Ind_q`, where :math:`Ind_p = [p_1, p_2]` and :math:`Ind_q = [q_1, q_2]`, We find the *best* model :math:`ARMA(p,q)` that fits to the data and estimates the corresponding coefficients. We proceed as follows: - the object *WhittleFactory* is created with either a specified order :math:`(p,q)` or a range :math:`Ind_p \times Ind_q`. By default, the Welch estimator (object *Welch*) is used with its default parameters. - for each order :math:`(p,q)`, the estimation of the parameters is done by maximizing the reduced equation of the Whittle likelihood function (\[lik2\]), thanks to the method *build* of the object *WhittleFactory*. This method applies to a time series or a process sample. If the user wants to get the quantified criteria :math:`AIC_c, AIC` and *BIC* of the model :math:`ARMA(p,q)`, he has to specify it by giving a *Point* of size 0 (*Point()*) as input parameter of the method *build*. - the output of the estimation is, in all the cases, one unique ARMA: the ARMA with the specified order or the optimal one with respect to the :math:`AIC_c` criterion. - in the case of a range :math:`Ind_p \times Ind_q`, the user can get all the estimated models thanks to the method *getHistory* of the object *WhittleFactory*. If the *build* has been parameterized by a *Point* of size 0, the user also has access to all the quantified criteria. The synthetic data is generated using the following 1-d ARMA process: .. math:: X_{0,t} + 0.4 X_{0,t-1} + 0.3 X_{0,t-2} + 0.2 X_{0,t-3} + 0.1 X_{0,t-4} = E_{0,t} + 0.4 E_{0,t-1} + 0.3 E_{0,t-2} with the noise :math:`E` defined as: .. math:: E \sim Triangular(-1, 0, 1) .. GENERATED FROM PYTHON SOURCE LINES 56-61 .. code-block:: Python import openturns as ot ot.RandomGenerator.SetSeed(0) ot.Log.Show(ot.Log.NONE) .. GENERATED FROM PYTHON SOURCE LINES 62-63 Create an arma process .. GENERATED FROM PYTHON SOURCE LINES 63-80 .. code-block:: Python tMin = 0.0 n = 1000 timeStep = 0.1 myTimeGrid = ot.RegularGrid(tMin, timeStep, n) myWhiteNoise = ot.WhiteNoise(ot.Triangular(-1.0, 0.0, 1.0), myTimeGrid) myARCoef = ot.ARMACoefficients([0.4, 0.3, 0.2, 0.1]) myMACoef = ot.ARMACoefficients([0.4, 0.3]) arma = ot.ARMA(myARCoef, myMACoef, myWhiteNoise) tseries = ot.TimeSeries(arma.getRealization()) # Create a sample of N time series from the process N = 100 sample = arma.getSample(N) .. GENERATED FROM PYTHON SOURCE LINES 81-82 CASE 1 : we specify a (p,q) order .. GENERATED FROM PYTHON SOURCE LINES 82-107 .. code-block:: Python # Specify the order (p,q) p = 4 q = 2 # Create the estimator factory = ot.WhittleFactory(p, q) print("Default spectral model factory = ", factory.getSpectralModelFactory()) # To set the spectral model factory # For example, set WelchFactory as SpectralModelFactory # with the Hann filtering window # The Welch estimator splits the time series in four blocs without overlap myFilteringWindow = ot.Hann() mySpectralFactory = ot.WelchFactory(myFilteringWindow, 4, 0) factory.setSpectralModelFactory(mySpectralFactory) print("New spectral model factory = ", factory.getSpectralModelFactory()) # Estimate the ARMA model from a time series # To get the quantified AICc, AIC and BIC criteria arma42, criterion = factory.buildWithCriteria(tseries) AICc, AIC, BIC = criterion[0:3] print("AICc=", AICc, "AIC=", AIC, "BIC=", BIC) arma42 .. rst-class:: sphx-glr-script-out .. code-block:: none Default spectral model factory = class=WelchFactory window = class=FilteringWindows implementation=class=Hamming blockNumber = 1 overlap = 0 New spectral model factory = class=WelchFactory window = class=FilteringWindows implementation=class=Hann blockNumber = 4 overlap = 0 AICc= 772.0387560411838 AIC= 771.0814910839188 BIC= 824.677883406151 .. raw:: html
class= ARMA timeGrid=class=RegularGrid name=Unnamed start=0 step=0.1 n=1000 coefficients AR=class=ARMACoefficients, shift=0, value=class=SquareMatrix dimension=1 implementation=class=MatrixImplementation name=Unnamed rows=1 columns=1 values=[-0.189305], shift=1, value=class=SquareMatrix dimension=1 implementation=class=MatrixImplementation name=Unnamed rows=1 columns=1 values=[0.424708], shift=2, value=class=SquareMatrix dimension=1 implementation=class=MatrixImplementation name=Unnamed rows=1 columns=1 values=[0.204215], shift=3, value=class=SquareMatrix dimension=1 implementation=class=MatrixImplementation name=Unnamed rows=1 columns=1 values=[0.0584364] coefficients MA=class=ARMACoefficients, shift=0, value=class=SquareMatrix dimension=1 implementation=class=MatrixImplementation name=Unnamed rows=1 columns=1 values=[-0.168415], shift=1, value=class=SquareMatrix dimension=1 implementation=class=MatrixImplementation name=Unnamed rows=1 columns=1 values=[0.452162] noiseDistribution= class=Normal name=Normal dimension=1 mean=class=Point name=Unnamed dimension=1 values=[0] sigma=class=Point name=Unnamed dimension=1 values=[0.406733] correlationMatrix=class=CorrelationMatrix dimension=1 implementation=class=MatrixImplementation name=Unnamed rows=1 columns=1 values=[1] state= class= ARMAState x= class=Sample name=Unnamed implementation=class=SampleImplementation name=Unnamed size=4 dimension=1 data=[[-0.0623912],[0.236768],[-0.107643],[0.10339]] epsilon= class=Sample name=Unnamed implementation=class=SampleImplementation name=Unnamed size=2 dimension=1 data=[[-0.26662],[0.0665265]]


.. GENERATED FROM PYTHON SOURCE LINES 108-109 CASE 2 : we specify a range of (p,q) orders .. GENERATED FROM PYTHON SOURCE LINES 109-125 .. code-block:: Python # Range for p pIndices = [1, 2, 4] # Range for q = [4,5,6] qIndices = [4, 5, 6] # Build a Whittle factory with default SpectralModelFactory (WelchFactory) # this time using ranges of order p and q factory_range = ot.WhittleFactory(pIndices, qIndices) # Estimate the arma model from a process sample arma_range, criterion = factory_range.buildWithCriteria(sample) AICc, AIC, BIC = criterion[0:3] print("AICc=", AICc, "AIC=", AIC, "BIC=", BIC) arma_range .. rst-class:: sphx-glr-script-out .. code-block:: none AICc= 4443.580404874942 AIC= 4443.35276259852 BIC= 4516.35727597643 .. raw:: html
class= ARMA timeGrid=class=RegularGrid name=Unnamed start=0 step=0.1 n=1000 coefficients AR=class=ARMACoefficients, shift=0, value=class=SquareMatrix dimension=1 implementation=class=MatrixImplementation name=Unnamed rows=1 columns=1 values=[0.419883], shift=1, value=class=SquareMatrix dimension=1 implementation=class=MatrixImplementation name=Unnamed rows=1 columns=1 values=[0.176036] coefficients MA=class=ARMACoefficients, shift=0, value=class=SquareMatrix dimension=1 implementation=class=MatrixImplementation name=Unnamed rows=1 columns=1 values=[0.422707], shift=1, value=class=SquareMatrix dimension=1 implementation=class=MatrixImplementation name=Unnamed rows=1 columns=1 values=[0.183076], shift=2, value=class=SquareMatrix dimension=1 implementation=class=MatrixImplementation name=Unnamed rows=1 columns=1 values=[-0.19104], shift=3, value=class=SquareMatrix dimension=1 implementation=class=MatrixImplementation name=Unnamed rows=1 columns=1 values=[-0.110271] noiseDistribution= class=Normal name=Normal dimension=1 mean=class=Point name=Unnamed dimension=1 values=[0] sigma=class=Point name=Unnamed dimension=1 values=[0.409622] correlationMatrix=class=CorrelationMatrix dimension=1 implementation=class=MatrixImplementation name=Unnamed rows=1 columns=1 values=[1] state= class= ARMAState x= class=Sample name=Unnamed implementation=class=SampleImplementation name=Unnamed size=2 dimension=1 data=[[-0.154588],[-0.0613563]] epsilon= class=Sample name=Unnamed implementation=class=SampleImplementation name=Unnamed size=4 dimension=1 data=[[-0.561741],[0.311484],[-0.0734148],[-0.179361]]


.. GENERATED FROM PYTHON SOURCE LINES 126-127 Results exploitation .. GENERATED FROM PYTHON SOURCE LINES 127-130 .. code-block:: Python # Get the white noise of the (best) estimated arma arma_range.getWhiteNoise() .. raw:: html
class=WhiteNoise distribution=class=Normal name=Normal dimension=1 mean=class=Point name=Unnamed dimension=1 values=[0] sigma=class=Point name=Unnamed dimension=1 values=[0.409622] correlationMatrix=class=CorrelationMatrix dimension=1 implementation=class=MatrixImplementation name=Unnamed rows=1 columns=1 values=[1]


.. _sphx_glr_download_auto_data_analysis_estimate_stochastic_processes_plot_estimate_arma.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_estimate_arma.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_estimate_arma.py `