# Inference for AR(p) Time Series

Consider a (stationary) autoregressive process, say of order 2,

for some white noise with variance . Here is a code to generate such a process,

> phi1=.25
> phi2=.7
> n=1000
> set.seed(1)
> e=rnorm(n)
> Z=rep(0,n)
> for(t in 3:n) Z[t]=phi1*Z[t-1]+phi2*Z[t-2]+e[t]
> Z=Z[800:1000]
> n=length(Z)
> plot(Z,type="l") Here, we have to estimate two sets of parameters: the autoregressive coefficients , and the variance of the innovation process . Several techniques can be used to estimate those parameters.

• using least square regression

A natural idea is to see here a regression model, since (if we consider a matrix formulation) Here we can run (conditional) ordinary least squares estimation,

> base=data.frame(Y=Z[3:n],X1=Z[2:(n-1)],X2=Z[1:(n-2)])
> regression=lm(Y~0+X1+X2,data=base)
> summary(regression)

Call:
lm(formula = Y ~ 0 + X1 + X2, data = base)

Residuals:
Min      1Q  Median      3Q     Max
-3.0268 -0.7063  0.1065  0.6925  3.2566

Coefficients:
Estimate Std. Error t value Pr(>|t|)
X1  0.23400    0.05463   4.283 2.88e-05 ***
X2  0.62863    0.05476  11.479  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.062 on 197 degrees of freedom
Multiple R-squared:  0.6349,	Adjusted R-squared:  0.6312
F-statistic: 171.3 on 2 and 197 DF,  p-value: < 2.2e-16

so we get the following estimators, for the autocorrelation coefficients, and the volatility of the noise

> regression$coefficients X1 X2 0.2339959 0.6286321 > summary(regression)$sigma
 1.061839
• using Yule-Walker equations

As we’ve seen in class, we can easily get the following equations for the autocovariance functions, which can also be written (again, using a matrix expression) So we just have to solve a simple linear system of equations. Note that if we divide by the variance , those equations can be written in terms of the autocorrelation functions The code is the following

> rho1=cor(Z[1:(n-1)],Z[2:n])
> rho2=cor(Z[1:(n-2)],Z[3:n])
> A=matrix(c(1,rho1,rho1,1),2,2)
> b=matrix(c(rho1,rho2),2,1)
> (PHI=solve(A,b))
[,1]
[1,] 0.2256270
[2,] 0.6315329

Now, we need to extract the estimated innovation process, from this set of parameters

> estWN=base$Y-(PHI*base$X1+PHI*base$X2) > sd(estWN)  1.058558 This estimator is probably not the best one (we can take into account that we’ve lost two degrees of freedom), but as a starting point, let us consider this one. An alternative could be to include the variance term in Yule-Walker equations, to get a three dimensional linear equation, It is not much more complicated to solve, actually, > gamma0=var(Z[1:n]) > gamma1=var(Z[1:(n-1)],Z[2:n]) > gamma2=var(Z[1:(n-2)],Z[3:n]) > A=matrix(c(gamma1,gamma0,gamma1,gamma2,gamma1,gamma0,1,0,0),3,3) > b=matrix(c(gamma0,gamma1,gamma2),3,1) > (PHISIGMA=solve(A,b)) [,1] [1,] 0.2283151 [2,] 0.6283431 [3,] 1.1335501 • using (conditional) likelihood estimators Finally, we can assume some distribution for the innovation process. The standard model is a Gaussian model, i.e. has a Gaussian distribution In that case, the conditional log likelihood (conditional since we set the first two observations here) is > CondLogLik=function(A,TS){ + phi1=A; phi2=A + sigma=A; L=0 + for(t in 3:length(TS)){ + L=L+dnorm(TS[t],mean=phi1*TS[t-1]+ + phi2*TS[t-2],sd=sigma,log=TRUE)} + return(-L)} Now, we can run standard optimization procedures, > LogL=function(A) CondLogLik(A,TS=Z) > optim(c(0,0,1),LogL)$par
 0.2339589 0.6285002 1.0565613

$value  293.3042$counts
106       NA

$convergence  0$message
NULL

It is also possible to consider a global maximum likelihood optimisation problem, since the variance matrix of vector has a know form.

• using (unconditional) likelihood estimators

The variance matrix of is , where autocovariances are not not know, be can easily be computed using a recursive relationship.

> library(mnormt)
> GlobalLogLik=function(A,TS){
+ n=length(TS)
+ phi1=A;  phi2=A
+ sigma=A
+ SIG=matrix(0,n,n)
+ rho=rep(0,n)
+ rho=1
+ rho=phi1/(1-phi2)
+ for(h in 3:n) rho[h]=phi1*rho[h-1]+phi2*rho[h-2]
+ for(i in 1:n){for(j in 1:n){
+ SIG[i,j]=rho[abs(i-j)+1]}}
+ gamma0=(1-phi2)*sigma^2/((1+phi2)*((1-phi2)^2-phi1^2))
+ SIG=gamma0*SIG
+ return(dmnorm(TS,rep(0,n),SIG,log=TRUE))}
> LogL=function(A) -GlobalLogLik(A,TS=Z)
> optim(c(.1,.1,1),LogL)
Error in chol.default(x, pivot = FALSE) :
Error in pd.solve(varcov, log.det = TRUE) :
x appears to be not positive definite

The problem is that there is a strong constraint on the pair to get a stationary process (we are not far away, here, from the border of the triangle, where the process become non stationary). To be more specific (this was mentioned in a previous post), we should have

i.e. in a standard matrix form

(we can add an additional constraint on the variance parameter, to insure that it will be positive). To run a contrained optimization routine, consider

> U=matrix(c(1,0,0,-1,0,1,0,-1,0,0,1,0),4,3)
> C=c(0,0,0,-.99999)
$par  0.2238892 0.6342850 1.0613388$value
 297.9202

$counts function gradient 108 NA$convergence
 0

$message NULL$outer.iterations
 2

$barrier.value  0.000189892 (here, to faster, we restrain the parameters so that they will be positive). • comparing those estimates Here, our five estimators are rather close. Let us run more samples to see more precisely how they behave. For the first parameter , we get and for the second one, , we have The bias we observe is probably coming from the fact that, with this numerical example, we are not far away from the non-stationary case (the sum of the true parameters should be less than 1, and it is 0.95). When we estimate the parameters, we force them to be inside the triangle, since those parameters can be estimated only if the process is stationary. Observe that the standard-deviation of the innovation process is here, well estimated, (with clearly some estimators that perform better than others). # Triangle for Parameters of AR(2) Stationary Processes We’ve seen yesterday conditions on so that the canonical process, , satisfying The condition is rather simple, since should be a triangular region. But the proof is a bit more tricky… Recall that we want to parametrize the region Since we have a true process, then . Our polynomial is here where ‘s are the roots – in – of . Consider now some kind of dual version of that polynomial, Having the roots of outside the unit circle is the same as having the roots of inside the unit circle. Obserse that we can write Roots of are then From this point, we should discuss a little bit, depending on the value of . • if Then there is one root, and only one. So we need to have or equivalently . • if Then we got roots in , and means, equivalently, that • if Then we have two (conjugate) roots in , and the square of norm of those roots is . Thus, . We get what was mention in the course: the canonical has a stationary solution if, and only if which is a triangular region, see # Visualizing Autoregressive Time Series In the MAT8181 graduate course on Time Series, we started discussing autoregressive models. Just to illustrate, here is some code to plot – causal – process, > graphar1=function(phi){ + nf <- layout(matrix(c(1,1,1,1,2,3,4,5), 2, 4, byrow=TRUE), respect=TRUE) + e=rnorm(n) + X=rep(0,n) + for(t in 2:n) X[t]=phi*X[t-1]+e[t] + plot(X[1:6000],type="l",ylab="") + abline(h=mean(X),lwd=2,col="red") + abline(h=mean(X)+2*sd(X),lty=2,col="red") + abline(h=mean(X)-2*sd(X),lty=2,col="red") + u=seq(-1,1,by=.001) + plot(0:1,0:1,col="white",xlab="",ylab="",axes=FALSE,ylim=c(-2,2),xlim=c(-2.5,2.5)) + polygon(c(u,rev(u)),c(sqrt(1-u^2),rev(-sqrt(1-u^2))),col="light yellow") + abline(v=0,col="grey") + abline(h=0,col="grey") + points(1/phi,0,pch=19,col="red",cex=1.3) + plot(0:1,0:1,col="white",xlab="",ylab="",axes=FALSE,ylim=c(-.2,.2),xlim=c(-1,1)) + axis(1) + points(phi,0,pch=19,col="red",cex=1.3) + acf(X,lwd=3,col="blue",main="",ylim=c(-1,1)) + pacf(X,lwd=3,col="blue",main="",ylim=c(-1,1),xlim=c(0,16))} e.g. > graphar1(.8) or > graphar1(-.7) (with, on the bottom, the root of the characteristic polynomial, the value of the parameter , the autocorrelation function and the partial autocorrelation function ). Of course, it is possible to do something similar with processes, > graphar2=function(phi1,phi2){ + nf <- layout(matrix(c(1,1,1,1,2,3,4,5), 2, 4, byrow=TRUE), respect=TRUE) + e=rnorm(n) + X=rep(0,n) + for(t in 3:n) X[t]=phi1*X[t-1]+phi2*X[t-2]+e[t] + plot(X[1:6000],type="l",ylab="") + abline(h=mean(X),lwd=2,col="red") + abline(h=mean(X)+2*sd(X),lty=2,col="red") + abline(h=mean(X)-2*sd(X),lty=2,col="red") + P=polyroot(c(1,-phi1,-phi2)) + u=seq(-1,1,by=.001) + plot(0:1,0:1,col="white",xlab="",ylab="",axes=FALSE,ylim=c(-2,2),xlim=c(-2.5,2.5)) + polygon(c(u,rev(u)),c(sqrt(1-u^2),rev(-sqrt(1-u^2))),col="light yellow") + abline(v=0,col="grey") + abline(h=0,col="grey") + points(P,pch=19,col="red",cex=1.3) + plot(0:1,0:1,col="white",xlab="",ylab="",axes=FALSE,xlim=c(-2.1,2.1),ylim=c(-1.2,1.2)) + polygon(c(-2,0,2,-2),c(-1,1,-1,-1),col="light green") + u=seq(-2,2,by=.001) + lines(u,-u^2/4) + abline(v=seq(-2,2,by=.2),col="grey",lty=2) + abline(h=seq(-1,1,by=.2),col="grey",lty=2) + segments(0,-1,0,1) + axis(1) + axis(2) + points(phi1,phi2,pch=19,col="red",cex=1.3) + acf(X,lwd=3,col="blue",main="",ylim=c(-1,1)) + pacf(X,lwd=3,col="blue",main="",ylim=c(-1,1),xlim=c(0,16))} For example, > graphar2(.65,.3) or > graphar2(-1.4,-.7) # Simulation de séries temporelles Un billet rapide pour reprendre le code tapé en cours, la semaine passée. Considérons un processus autorégressif d’ordre 1, où est un bruit blanc, stationnaire, i.e. appartient à l’intervalle . Le code pour simuler un tel processus est n=1000 bruit=rnorm(n) phi1= .85 X=rep(NA,n) X=0 for(t in 2:n){X[t]=phi1*X[t-1]+bruit[t]} plot(acf(X),lwd=5,col='blue') plot(pacf(X),lwd=5,col='blue')  ou avec un autocorrélation au premier ordre négative, phi1= -0.7  On peut aussi regarder un processus autorégressif au second ordre, sur la figure ci-dessous (avec en haut à gauche le triangle de stationnarité du couple de paramètres). phi1= 0.3 phi2= 0.5 X=rep(NA,n) X[1:2]=0 for(t in 3:n){ X[t]=phi1*X[t-1]+phi2*X[t-2]+bruit[t]} Histoire de changer un peu, on peut regarder un processus moyenne mobile au premier ordre, où est un paramètre dans . theta1= .8 X=rep(NA,n) X=0 for(t in 2:n){ X[t]=bruit[t]+theta1*bruit[t-1]} ou une moyenne mobile du second ordre, theta1= -.6 theta2= .5 X=rep(NA,n) X[1:2]=0 for(t in 3:n){ X[t]=bruit[t]+theta1*bruit[t-1]+ theta2*bruit[t-2]} # Inference and autoregressive processes Consider a (stationary) autoregressive process, say of order 2, for some white noise with variance . Here is a code to generate such a process, > phi1=.5 > phi2=-.4 > sigma=1.5 > set.seed(1) > n=240 > WN=rnorm(n,sd=sigma) > Z=rep(NA,n) > Z[1:2]=rnorm(2,0,1) > for(t in 3:n){Z[t]=phi1*Z[t-1]+phi2*Z[t-2]+WN[t]} Here, we have to estimate two sets of parameters: the autoregressive coefficients , and the variance of the innovation process . There are (at least) three techniques to estimate those parameters. • using least square regression A natural idea is to see here a regression model, and thus, if we consider a matrix formulation, Here we can run (conditional) ordinary least squares estimation, > base=data.frame(Y=Z[3:n],X1=Z[2:(n-1)],X2=Z[1:(n-2)]) > regression=lm(Y~0+X1+X2,data=base) > summary(regression) Call: lm(formula = Y ~ 0 + X1 + X2, data = base) Residuals: Min 1Q Median 3Q Max -4.3491 -0.8890 -0.0762 0.9601 3.6105 Coefficients: Estimate Std. Error t value Pr(>|t|) X1 0.45107 0.05924 7.615 6.34e-13 *** X2 -0.41454 0.05924 -6.998 2.67e-11 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.449 on 236 degrees of freedom Multiple R-squared: 0.2561, Adjusted R-squared: 0.2497 F-statistic: 40.61 on 2 and 236 DF, p-value: 6.949e-16 > regression$coefficients
X1         X2
0.4510703 -0.4145365
> summary(regression)$sigma  1.449276 • using Yule-Walker equations As we’ve seen in class, we can easily get the following equations for the autocovariance functions, which can also be written So we just have to solve a simple linear system of equations. Note that if we divide by the variance , those equations can be written in terms of the autocorrelation functions The code is the following > rho1=cor(Z[1:(n-1)],Z[2:n]) > rho2=cor(Z[1:(n-2)],Z[3:n]) > A=matrix(c(1,rho1,rho1,1),2,2) > b=matrix(c(rho1,rho2),2,1) > (PHI=solve(A,b)) [,1] [1,] 0.4517579 [2,] -0.4155920 Now, we need to extract the estimated innovation process, from this set of parameters (note that it could be possible to include the variance term in Yule-Walker equations, to get a three dimensional linear equation) > estWN=base$Y-(PHI*base$X1+PHI*base$X2)
> sd(estWN)
 1.445706

This estimator is probably not the best one (we can take into account that we’ve lost two degrees of freedom), but as a starting point, let us consider this one.

• using (conditional) likelihood estimators

Finally, we can assume some distribution for the innovation process. Thestandard model is a Gaussian model, i.e. In that case, the conditional log likelihood (conditional since we set the first two observations here) is

> CondLogLik=function(A,TS){
+ phi1=A;  phi2=A
+ sigma=A	; L=0
+ for(t in 3:length(TS)){
+ L=L+dnorm(TS[t],mean=phi1*TS[t-1]+
+ phi2*TS[t-2],sd=sigma,log=TRUE)}
+ return(-L)}

Now, we can run standard optimization procedures,

> LogL=function(A) CondLogLik(A,TS=Z)
> optim(c(0,0,1),LogL)
$par  0.4509685 -0.4144938 1.4430930$value
 425.0164

$counts function gradient 88 NA$convergence
 0

\$message
NULL

Here, our three estimators are rather close. Actually, if we generate 1,000 time series (of size 240), those are the Box-plots of our three estimators, for the first order autoregressive coefficient for the second one, and finally for the standard deviation of the innovation process All those estimators behave nicely, and are rather close. Note that they all might be biased, but they are consistent (see Davidson and MacKinnon for instance, in their book, for more details).