Estimate a scalar ARMA process¶

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 $(p,q)$ , the library fits a model $ARMA(p,q)$ to the data by estimating the coefficients $(a_1, \dots, a_p), (b_1, \dots, b_q)$ and the variance $\sigma$ of the white noise.

If the User specifies a range of orders $(p,q) \in Ind_p \times Ind_q$ , where $Ind_p = [p_1, p_2]$ and $Ind_q = [q_1, q_2]$ , We find the best model $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 $(p,q)$ or a range $Ind_p \times Ind_q$ . By default, the Welch estimator (object Welch) is used with its default parameters.
for each order $(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 $AIC_c, AIC$ and BIC of the model $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 $AIC_c$ criterion.
in the case of a range $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:

$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 $E$ defined as:

$E \sim Triangular(-1, 0, 1)$

import openturns as ot

ot.RandomGenerator.SetSeed(0)
ot.Log.Show(ot.Log.NONE)

Create an arma process

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)

CASE 1 : we specify a (p,q) order

# 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

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

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]]

CASE 2 : we specify a range of (p,q) orders

# 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

AICc= 4443.580404874942 AIC= 4443.35276259852 BIC= 4516.35727597643

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]]

Results exploitation

# Get the white noise of the (best) estimated arma
arma_range.getWhiteNoise()

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]

OpenTURNS

An Open source initiative for the Treatment of Uncertainties, Risks'N Statistics

Previous topic

Next topic

This Page

Estimate a scalar ARMA process¶