Après les exposés des dernières séances, l’examen du cours MAT8181, **Séries Chronologiques** avait lieu ce matin (et devrait finir dans quelques minutes, avec un peu de temps supplémentaire pour certain, compte tenu de la panne de métro qu’on a eu la joie de subir). L’énoncé est en ligne, et j’ai aussi écrit quelques éléments de correction. En cas de désaccord (mineur ou majeur) avec mes réponses, merci de me le faire savoir rapidement !

# Tag Archives: MAT8181

# Inference for ARCH processes

Consider some ARCH() process, say ARCH(),

where

with a Gaussian (strong) white noise .

> n=500 > a1=0.8 > a2=0.0 > w= 0.2 > set.seed(1) > eta=rnorm(n) > epsilon=rnorm(n) > sigma2=rep(w,n) > for(t in 3:n){ + sigma2[t]=w+a1*epsilon[t-1]^2+a2*epsilon[t-2]^2 + epsilon[t]=eta[t]*sqrt(sigma2[t]) + } > par(mfrow=c(1,1)) > plot(epsilon,type="l",ylim=c(min(epsilon)-.5,max(epsilon))) > lines(min(epsilon)-1+sqrt(sigma2),col="red")

(the red line is the conditional variance process).

> par(mfrow=c(1,2)) > acf(epsilon,lag=50,lwd=2) > acf(epsilon^2,lag=50,lwd=2)

We did mention in class that if a ARCH(), then is an AR() process. So a first idea is to consider a regression, as we did for Gaussian AR()

> db=data.frame(Y=epsilon[2:n]^2,X1=epsilon[1:(n-1)]^2) > summary(lm(Y~X1,data=db)) Call: lm(formula = Y ~ X1, data = db) Residuals: Min 1Q Median 3Q Max -2.4538 -0.3618 -0.2626 0.0935 9.3667 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.34963 0.04342 8.052 6.08e-15 *** X1 0.31123 0.04262 7.303 1.13e-12 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.8413 on 497 degrees of freedom Multiple R-squared: 0.0969, Adjusted R-squared: 0.09508 F-statistic: 53.33 on 1 and 497 DF, p-value: 1.129e-12

There is some significant autocorrelation here. But since our vectors cannot be considered as Gaussian, using least squares is perhaps not the best strategy. Actually, if our series is not Gaussian, it is still conditionally Gaussian, since we assumed that is a Gaussian (strong) white noise,

The likelihood is then

and the log-likelihood is

And a natural idea is to define

The code is simply

> X=epsilon > loglik=function(param){ + w=exp(param[1]) + a1=exp(param[2]) + s2=rep(w,n) + for(t in 2:length(X)){s2[t]=w+a1*X[t-1]^2} + logL=-.5*sum(log(s2))-.5*sum(X^2/s2) + return(-logL) + } > OPT=optim(par= + coefficients(lm(Y~X1,data=db)),fn=loglik) > exp(OPT$par) (Intercept) X1 0.2482241 0.5858578

(since the parameters have to be positive, we assume here that they can be written as the exponential of some real values). Observe that those values are closer to the one used to generate our time series.

If we use R functions to estimate those parameters, we get

> library(tseries) > summary(garch(epsilon,c(0,1))) ... Call: garch(x = epsilon, order = c(0, 1)) Model: GARCH(0,1) Residuals: Min 1Q Median 3Q Max -2.87023 -0.60836 -0.03426 0.66648 3.48443 Coefficient(s): Estimate Std. Error t value Pr(>|t|) a0 0.24959 0.02470 10.104 < 2e-16 *** a1 0.58306 0.09737 5.988 2.13e-09 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

so that the confidence interval for is

> summary(garch(epsilon,c(0,1)))$coef[2,1]+ + c(-1.96,1.96)*summary(garch(epsilon,c(0,1)))$coef[2,2] [1] 0.3922088 0.7739088

Actually, since our main interest is this parameter, it is possible to use profile likelihood techniques,

> proflik=function(a){ + loglik=function(w){ + s2=rep(w,n) + for(t in 2:length(X)){s2[t]=w+a*X[t-1]^2} + logL=-.5*sum(log(s2))-.5*sum(X^2/s2) + return(-logL)} + return(-optim(par=.3,fn=loglik)$value)} > A=seq(0,2,by=.05) > P=Vectorize(proflik)(A) > par(mfrow=c(1,1)) > plot(A,P,type="l") > OPT=optimize(function(x) -proflik(x), interval=c(0,2)) > t=-OPT$objective-qchisq(.95,df=1) > abline(h=t,col="red") > ainf=uniroot(function(x) proflik(x)-t,c(0,OPT$minimum))$root > asup=uniroot(function(x) proflik(x)-t,c(OPT$minimum,2))$root > abline(v=ainf,lty=2) > abline(v=asup,lty=2)

Of course, all those techniques can be extended to higher order ARCH processes. For instance, if we assume that we have a ARCH() time series

where now

with a Gaussian (strong) white noise . The log-likelihood is still

and we can define

The code above can be changed, to take into account this additional component,

> db=data.frame(Y=epsilon[3:n]^2, + X1=epsilon[2:(n-1)]^2, + X2=epsilon[1:(n-2)]^2) > X=epsilon > loglik=function(param){ + w=exp(param[1]) + a1=exp(param[2]) + a2=exp(param[3]) + s2=rep(w,n) + for(t in 3:length(X)){s2[t]=w+a1*X[t-1]^2+a2*X[t-2]^2} + logL=-.5*sum(log(s2))-.5*sum(X^2/s2) + return(-logL) + } > OPT=optim(par= + coefficients(lm(Y~X1+X2,data=db)),fn=loglik) > exp(OPT$par) (Intercept) X1 X2 0.22710526 0.59475474 0.04741294

We can also consider some Generalized ARCH process, e.g. a GARCH(,),

where now

Again, maximum likelihood techniques can be used. Actually, we can also code Fisher-Scoring algorithm, since (in a very general context)

with here . Using a standard gradient descent algorithm, we get the following estimate for our GARCH process,

> X=epsilon > theta=c(.2,.2,.2) > G=rep(1,3) > n=length(X) > j=1 > while(sum(G^2)>1e-12){ + s2=rep(theta[1],n) + for (i in 2:n){s2[i]=theta[1]+theta[2]*X[(i-1)]^2+theta[3]*s2[(i-1)]} + z=(X^2-s2)/s2^2 + V=cbind(z[2:n],z[2:n]*X[1:(n-1)]^2,z[2:n]*s2[1:(n-1)]) + H=(t(V)%*%V) + G=apply(V,2,sum) + theta=theta+solve(H)%*%G + j=j+1} > as.numeric(theta) [1] 0.20372918 0.59183911 0.08936159

The interesting point, here, is that we also derive the (asymptotic) variance

> (stdev=sqrt(diag(solve(H)))) [1] 0.01849067 0.04950477 0.02937233

# Independence and correlation

A short post to get back on a property I gave briefly in the MAT8595 class in January, and again in the MAT8181 class this week (to illustrate the distinction between weak and strong white noises). Recall that (real-valued) random variables and are independent if for all , Another characterization, for integrable variable is that for all , which can be written, if variables are square integrable The idea to prove this characterization is to observe that if and are independent can be written Using a standard argument in integration theory, equality is valid for step functions (not only indicators), and then to positive measurable functions, and finally to integrable functions. Proving this result is not that difficult. Observe that Rényi (1959) – inspired by Gebelein (1947) – followed by Sarmanov (1958) introduced the concept of *maximal correlation*, that can be related to this result, where the maximum is taken over all functions and such that the correlation exist. Actually, it is possible to consider only transformations such that and (and similarly for , the idea is that we simple center and scale, which does not impact the correlation.Thus, and are independent if and only if Algorithm to estimate that coefficient are interesting. The problem can be written, equivalently And if the minimization is considered over , assuming that is fixed, then the optimal transformation is And similarly for . So using an iterative algorithm, it is possible to get and (see Breiman and Friedman (1985) for more details). Actually, those functions appear in nonlinear canonical analysis. As mentioned in Lancaster (1957), for a Gaussian random vector and in that case and are affine functions. This can be related to Hermite’s polynomial and to the expansion of the bivariate Gaussian density. I still hope that someone will go further for the project in the MAT8181 course.

# Seasonal Unit Roots

As discussed in the MAT8181 course, there are – at least – two kinds of non-stationary time series: those with a trend, and those with a unit-root (they will be called *integrated*). Unit root tests cannot be used to assess whether a time series is stationary, or not. They can only detect integrated time series. And the same holds for *seasonal unit root*.

In a previous post, we’ve seen that it was difficult to model hourly temperature, since most test do not reject unit roots. Consider here the average *monthly* temperature, still in Montréal, QC.

> montreal=read.table("http://freakonometrics.free.fr/temp-montreal-monthly.txt") > M=as.matrix(montreal[,2:13]) > X=as.numeric(t(M)) > tsm=ts(X,start=1948,freq=12) > plot(tsm)

For those who don’t know Montréal, Winter and Summer are very different. We can visualize monthly differences using

> month=rep(1:12,length(tsm)/12) > plot(month,as.numeric(tsm)) > lines(1:12,apply(M,2,mean),col="red",type="b",pch=19)

or, if install the uroot package, removed from the CRAN repository, we can use

> library(uroot) > bbplot(tsm)

or

> bb3D(tsm) Loading required package: tcltk

It looks like our time series is cyclic, because of the yearly seasonal pattern. The autocorrelation function is here

> acf(tsm,lag=120)

Again, this cycle can be visualized using

> persp(1948:2013,1:12,M,theta=-50,col="yellow",shade=TRUE, + xlab="Year",ylab="Month",zlab="Temperature",ticktype="detailed")

Now, the question is *is there a seasonal unit root *? This would mean that our model should be

If we forget about the autoregressive and the moving average component, we can estimate

If there is a *seasonal unit root* then should be close to 1. Somehow.

> arima(tsm,order=c(0,0,0),seasonal=list(order=c(1,0,0),period=12)) Call: arima(x = tsm, order = c(0, 0, 0), seasonal = list(order = c(1, 0, 0), period = 12)) Coefficients: sar1 intercept 0.9702 6.4566 s.e. 0.0071 2.1515

It is not far away from 1. Actually, it cannot be too close to 1. If it was, then we would get an error message…

To illustrate some interesting models, let us consider also *quarterly* temperatures,

> N=cbind(apply(montreal[,2:4],1,sum),apply(montreal[,5:7],1,sum),apply(montreal[,8:10],1,sum),apply(montreal[,11:13],1,sum)) > X=as.numeric(t(N)) > tsq=ts(X,start=1948,freq=4) > persp(1948:2013,1:4,N,theta=-50,col="yellow",shade=TRUE, + xlab="Year",ylab="Quarter",zlab="Temperature",ticktype="detailed")

(again, the aim is just to be able to write down some equations, if necessary)

Why not consider a model on the *quarterly *temperature? Something like

i.e.

where is some matrix . This model can easily be estimated,

> library(vars) > df=data.frame(N) > names(df)=paste("y",1:4,sep="") > model=VAR(df) > model VAR Estimation Results: ======================= Estimated coefficients for equation y1: ======================================= Call: y1 = y1.l1 + y2.l1 + y3.l1 + y4.l1 + const y1.l1 y2.l1 y3.l1 y4.l1 const -0.13943065 0.21451118 0.08921237 0.30362065 -34.74793931 Estimated coefficients for equation y2: ======================================= Call: y2 = y1.l1 + y2.l1 + y3.l1 + y4.l1 + const y1.l1 y2.l1 y3.l1 y4.l1 const 0.02520938 0.05288958 -0.13277377 0.05134148 40.68955266 Estimated coefficients for equation y3: ======================================= Call: y3 = y1.l1 + y2.l1 + y3.l1 + y4.l1 + const y1.l1 y2.l1 y3.l1 y4.l1 const 0.07740824 -0.21142726 0.11180518 0.12963931 56.81087283 Estimated coefficients for equation y4: ======================================= Call: y4 = y1.l1 + y2.l1 + y3.l1 + y4.l1 + const y1.l1 y2.l1 y3.l1 y4.l1 const 0.18842863 -0.31964579 0.25099508 -0.04452577 5.73228873

and matrix is here

> A=rbind( + coefficients(model$varresult$y1)[1:4], + coefficients(model$varresult$y2)[1:4], + coefficients(model$varresult$y3)[1:4], + coefficients(model$varresult$y4)[1:4]) > A y1.l1 y2.l1 y3.l1 y4.l1 [1,] -0.13943065 0.21451118 0.08921237 0.30362065 [2,] 0.02520938 0.05288958 -0.13277377 0.05134148 [3,] 0.07740824 -0.21142726 0.11180518 0.12963931 [4,] 0.18842863 -0.31964579 0.25099508 -0.04452577

Since stationary of this multiple time series is closely related to eigenvalues of this matrix, let us look at them,

> eigen(A)$values [1] 0.35834830 -0.32824657 -0.14042175 0.09105836 > Mod(eigen(A)$values) [1] 0.35834830 0.32824657 0.14042175 0.09105836

So it looks like there is no stationarity issue, here. A restricted model is the *periodic autoregressive model*, so called model, discussed by Paap and Franses

where

and

Keep in mind that this is a model, since

This model can be estimated using a specific package (one can also look at the vignette, to get a better understanding of the syntax)

> library(partsm) > detcomp <- list(regular=c(0,0,0), seasonal=c(1,0), regvar=0) > model=fit.ar.par(wts=tsq, detcomp=detcomp, type="PAR", p=1) > PAR.MVrepr(model) ---- Multivariate representation of a PAR model. Phi0: 1.000 0.000 0.000 0 -0.242 1.000 0.000 0 0.000 -0.261 1.000 0 0.000 0.000 -0.492 1 Phi1: 0 0 0 0.314 0 0 0 0.000 0 0 0 0.000 0 0 0 0.000 Eigen values of Gamma = Phi0^{-1} %*% Phi1: 0.01 0 0 0 Time varing accumulation of shocks: 0.010 0.040 0.155 0.314 0.002 0.010 0.037 0.076 0.001 0.003 0.010 0.020 0.000 0.001 0.005 0.010

Here, the characteristic equation is

so there is a (seasonal) unit root if

Which is clearly not the case, here. It is possible to perform Canova-Hansen test,

> CH.test(tsq) ------ - ------ ---- Canova & Hansen test ------ - ------ ---- Null hypothesis: Stationarity. Alternative hypothesis: Unit root. Frequency of the tested cycles: pi/2 , pi , L-statistic: 1.122 Lag truncation parameter: 5 Critical values: 0.10 0.05 0.025 0.01 0.846 1.01 1.16 1.35

The idea is that polynomial has four root, in ,

since

If we get back to monthly data, has twelve roots,

each of them having different interpretations.

Here we can have 1 cycle per year (on 12 months), 2 cycles per year (on 6 months), 3 cycles per year (on 4 months), 4 cycles per year (on 3 months), even 6 cycles per year (on 2 months). This will depend on the argument of the root, with respectively

The output of the test is here

> CH.test(tsm) ------ - ------ ---- Canova & Hansen test ------ - ------ ---- Null hypothesis: Stationarity. Alternative hypothesis: Unit root. Frequency of the tested cycles: pi/6 , pi/3 , pi/2 , 2pi/3 , 5pi/6 , pi , L-statistic: 1.964 Lag truncation parameter: 20 Critical values: 0.10 0.05 0.025 0.01 2.49 2.75 2.99 3.27

And it looks like we reject the assumption of a *seasonal *unit root. I can even mention the following testing procedure

> library(forecast) > nsdiffs(tsm, test="ch") [1] 0

where the ouput: “1” means that there is a seasonal unit root and “0” that there is no seasonal unit root. Simple to read, isn’t it? If we consider the periodic autoregressive model on the monthly data, the output is

> model=fit.ar.par(wts=tsm, detcomp=detcomp, type="PAR", p=1) > model ---- PAR model of order 1 . y_t = alpha_{1,s}*y_{t-1} + alpha_{2,s}*y_{t-2} + ... + alpha_{p,s}*y_{t-p} + coeffs*detcomp + epsilon_t, for s=1,2,...,12 ---- Autoregressive coefficients. s=1 s=2 s=3 s=4 s=5 s=6 s=7 s=8 s=9 s=10 s=11 s=12 alpha_1s 0.15 0.05 0.07 0.33 0.11 0 0.3 0.38 0.31 0.19 0.15 0.37

So, whatever the test, we always reject the assumption that there is a seasonal unit root. Which does not mean that we can not have a strong cycle! Actually, the series is almost periodic. But there is no unit root! So all of this makes sense (I hardly believe that there might be unit root – seasonal, or not – in temperatures).

Just to make sure that we get it right, consider two times series, inspired from the previous one. The first one is a periodic sequence (with a very very small noise, just to avoid problem of non-definite matrices) and the second one is clearly integrated.

> Xp1=Xp2=as.numeric(t(M)) > for(t in 13:length(M)){ + Xp1[t]=Xp1[t-12] + Xp2[t]=Xp2[t-12]+rnorm(1,0,2) + } > Xp1=Xp1+rnorm(length(Xp1),0,.02) > tsp1=ts(Xp1,start=1948,freq=12) > tsp2=ts(Xp2,start=1948,freq=12) > par(mfrow=c(2,1)) > plot(tsp1) > plot(tsp2)

see also

> par(mfrow=c(1,2)) > bb3D(tsp1) > bb3D(tsp2)

If we quickly look at those series, I would say that the first one has no unit root – even if it is not stationary, but it is because the series is periodic – while there is (are ?) unit root(s) for the second one. If we look at Canova-Hansen test, we get

> CH.test(tsp1) ------ - ------ ---- Canova & Hansen test ------ - ------ ---- Null hypothesis: Stationarity. Alternative hypothesis: Unit root. Frequency of the tested cycles: pi/6 , pi/3 , pi/2 , 2pi/3 , 5pi/6 , pi , L-statistic: 2.234 Lag truncation parameter: 20 Critical values: 0.10 0.05 0.025 0.01 2.49 2.75 2.99 3.27 > CH.test(tsp2) ------ - ------ ---- Canova & Hansen test ------ - ------ ---- Null hypothesis: Stationarity. Alternative hypothesis: Unit root. Frequency of the tested cycles: pi/6 , pi/3 , pi/2 , 2pi/3 , 5pi/6 , pi , L-statistic: 5.448 Lag truncation parameter: 20 Critical values: 0.10 0.05 0.025 0.01 2.49 2.75 2.99 3.27

I know that this package has been removed, so maybe I should not use it. Consider instead

> nsdiffs(tsp1, 12,test="ch") [1] 0 > nsdiffs(tsp2, 12,test="ch") [1] 1

Here we have the same conclusion. The first one does not have a unit root, but the second one has. But be careful: with Osborn-Chui-Smith-Birchenhall test,

> nsdiffs(tsp1, 12,test="ocsb") [1] 1 > nsdiffs(tsp2, 12,test="ocsb") [1] 1

we have the feeling that there is also a unit root in our cyclic series…

So here, on a short-frequency, we do reject the assumption of a unit root – even a seasonal one – on our temperature series. We still have our high-frequency problem to deal with, some day (but I don’t think I’ll have enough time to introduce long range dependence this session, unfortunately).

# Linear ‘Prediction’ for AR Time Series

In the exercises for the MAT8181 course, there are two Exercises (16 and 17) about prediction and extrapolation based on MA(1) and AR(1) time series. But before discussing those exercises (I had some request for hints), I wanted to recall the definition of the linear prediction,

where

As discussed previously on this blog, we consider here on projection not on (this would be the conditional expectancy) but on the linear subset.

The goal of Exercise 2 was to establish an important result, in the context of Gaussian random vectors. If is a (multivariate) Gaussian vector, , then

where is the vector .

Keep those results in mind, and let us look at Exercise 17, for instance. Here, is an AR(1) process, with innovation ,

One observation (say ) is missing. We have here 3 questions:

- what is the best linear prediction of given and
- what is the best linear prediction of given and
- what is the best linear prediction of given and

**Case 1.** Here, we have to compute

Since we have an AR(1) process, , and . Thus, from the relationship above

which can be written

i.e. . Which makes sense actually: the AR(1) process is Markovian of order one, so

And we seen in class that for an AR(1) process

So far, so good.

**Case 2.** Here, we have to compute

Since we have an AR(1) process, , and . Thus, from the relationship above

i.e. .

**Case 3.** Finally, we have to compute

One more time, since we have an AR(1) process, , and . So here, the relationship above becomes

Here, we can write

i.e.

So finally, what we got here is

and

The mean squared errors for each of those estimates are obtained computing

I guess I should probably stop here… that’s a detailed hint actually.

# Seasonal, or periodic, time series

Monday, in our MAT8181 class, we’ve discussed seasonal unit roots from a practical perspective (the theory will be briefly mentioned in a few weeks, once we’ve seen multivariate models). Consider some time series , for instance traffic on French roads,

> autoroute=read.table( + "http://freakonometrics.blog.free.fr/public/data/autoroute.csv", + header=TRUE,sep=";") > X=autoroute$a100 > T=1:length(X) > plot(T,X,type="l",xlim=c(0,120)) > reg=lm(X~T) > abline(reg,col="red")

As discussed in a previous post, if there is a trend, we should remove it, and work on the residual

> Y=residuals(reg) > acf(Y,lag=36,lwd=3)

We can observe that there is some seasonality, here. A first strategy might be to assume that there is a seasonal unit root, so we consider , and we try to find some ARMA process. Consider the empirical autocorrelation function of that time series,

> Z=diff(Y,12) > acf(Z,lag=36,lwd=3)

or the partial autocorrelation function

> pacf(Z,lag=36,lwd=3)

The first graph might suggest a MA(1) structure, while the second graph might suggest an AR(1) time series. Let us try both.

> model1=arima(Z,order=c(0,0,1)) > model1 Call: arima(x = Z, order = c(0, 0, 1)) Coefficients: ma1 intercept -0.2367 -583.7761 s.e. 0.0916 254.8805 sigma^2 estimated as 8071255: log likelihood = -684.1, aic = 1374.2 > E1=residuals(model1) > acf(E1,lag=36,lwd=3)

which can be considered as a white noise (if you’re not convinced, try either Box-Pierce or Ljung-Box test). Similarly,

> model2=arima(Z,order=c(1,0,0)) > model2 Call: arima(x = Z, order = c(1, 0, 0)) Coefficients: ar1 intercept -0.3214 -583.0943 s.e. 0.1112 248.8735 sigma^2 estimated as 7842043: log likelihood = -683.07, aic = 1372.15 > E2=residuals(model2) > acf(E2,lag=36,lwd=3)

which can be also considered as a white noise. So what we have, so far is

for some white noise . This suggest the following SARIMA structure on ,

> model2b=arima(Y,order=c(1,0,0), + seasonal = list(order = c(0, 1, 0), + period=12)) > model2b Call: arima(x = Y, order = c(1, 0, 0), seasonal = list(order = c(0, 1, 0), period = 12)) Coefficients: ar1 -0.2715 s.e. 0.1130 sigma^2 estimated as 8412999: log likelihood = -685.62, aic = 1375.25

So far, so good. Now, what if we consider that we do not have a seasonal unit root, but simply a large autoregressive coefficient in some AR structure. Let us try something like

where a natural guess is that this coefficient should – probably – be close to one. Let us try this one,

> model3c=arima(Y,order=c(1,0,0), + seasonal = list(order = c(1, 0, 0), + period = 12)) > model3c Call: arima(x = Y, order = c(1, 0, 0), seasonal = list(order = c(1, 0, 0), period = 12)) Coefficients: ar1 sar1 intercept -0.1629 0.9741 -684.9455 s.e. 0.1170 0.0115 3064.4040 sigma^2 estimated as 8406080: log likelihood = -816.11, aic = 1640.21

which is comparable with what we got previously (somehow), so we might assume that this model can be considered as an interesting one. We will discuss further the fact that the first coefficient might be considered as non-significant.

What is the difference from those two models? With a short term horizon, the two models are comparable. Clearly

> library(forecast) > previ=function(model,h=36,b=40000){ + prev=forecast(model,h) + T=1:85 + Tfutur=86:(85+h) + plot(T,Y,type="l",xlim=c(0,85+h),ylim=c(-b,b)) + polygon(c(Tfutur,rev(Tfutur)),c(prev$lower[,2],rev(prev$upper[,2])),col="orange",border=NA) + polygon(c(Tfutur,rev(Tfutur)),c(prev$lower[,1],rev(prev$upper[,1])),col="yellow",border=NA) + lines(prev$mean,col="blue") + lines(Tfutur,prev$lower[,2],col="red") + lines(Tfutur,prev$upper[,2],col="red") + }

Now, on a (very) long term perspective, the models are quite different: one is stationnary, so the forecast will tend to the average value (here 0, since the trend was removed), while the other one is (seasonaly) integrated, so the confidence interval will increase. For the non stationry, we get

> previ(model2b,600,b=60000)

and for the stationary one

> previ(model3c,600,b=60000)

But as mentioned in the introduction of this course, forecasts with those models are relevent only for short-term horizon (say not too large). And in that case, the prediction is almost the same here,

> previ(model2b,36,b=60000)

> previ(model3c,36,b=60000)

Now, if we come back on our second model, we did mention previously that the autoregressive coefficient might be considered as non-significant. What if we remove it?

> model3d=arima(Y,order=c(0,0,0), + seasonal = list(order = c(1, 0, 0), + period = 12)) > (model3d) Call: arima(x = Y, order = c(0, 0, 0), seasonal = list(order = c(1, 0, 0), period = 12)) Coefficients: sar1 intercept 0.9662 -696.5661 s.e. 0.0134 3182.3017 sigma^2 estimated as 8918630: log likelihood = -817.03, aic = 1640.07

If we look at a (short-term) forecast, we get

> previ(model3d,36,b=32000)

Do you see any difference? To be honest, I don’t… If we look at the figures, we get

> cbind(forecast(model2b,12)$mean,forecast(model3c,12)$mean,forecast(model3d,12)$mean) Time Series: Start = 86 End = 97 Frequency = 1 1 -4908.4920 -5092.8999 -5520.8780 2 -10012.7837 -9640.8103 -9493.0339 3 -3880.2202 -3841.1960 -3828.2611 4 -18102.5211 -17638.4086 -17499.1828 5 -20602.7346 -20090.9117 -19934.1066 6 -10463.2212 -10209.0139 -10132.0439 7 2458.1538 2376.4897 2351.2377 8 -1680.3342 -1654.4844 -1647.0057 9 876.6837 836.2342 823.4934 10 18046.5642 17561.6520 17413.1463 11 21531.4820 20956.3451 20780.2836 12 -3217.6103 -3152.0446 -3132.4112

Figures are different, but not significantly (keep in mind the size of the confidence interval). This might explain why, in R, when we ask for an autoregressive process or order , then we get a model with parameters to estimate, and even if some are not significant, we usually keep them for the forecast. Most of the time, from forecasting point of view, it’s no big deal.

# Filtering a Stationary Time Series

In the first part of the MAT8181 course, on linear (univariate) time series, I forgot to mention an important theorem. Let be a stationary time series, and a sequence of real numbers such that

then the time series defined as

is a stationary time series. Further, one can get easily that

This result can be used, if necessary in the exercises (that might save some time actually). I did not include this property in my notes because it is a bit technical to establish that this sum exists, and that the time series is stationary. It is rather simple with the spectral density (since where stands for the filter generating function), but I did not mention the spectral density since it requires some knowledge on Fourier analysis…

# Identification of ARMA processes

Last week (in the MAT8181 course) in order to identify the orders of an ARMA process, we’ve seen the eacf method, and I mentioned the scan method, introduced in Tsay and Tiao (1985). The code below – to produce the output of the scan procedure – has been adapted from an old code by Steve Chen (where I included a visualization of the *p*-values, with the following colors)

The procedure was described in the course, last Thursday,

arma.scan=function(z,ar.max=15,ma.max=15,alpha=0.01) { ym=function(z,t,m){return(z[t:(t-m)])} n=length(z) z=z - mean(z) cmax=ma.max + 1 rmax=ar.max + 1 corref=matrix(0,nrow=rmax,ncol=cmax) cmj.table=matrix(0,nrow=rmax,ncol=cmax) pv=matrix(0,nrow=rmax,ncol=cmax) mark=matrix(rep("X",(rmax)*(cmax)),nrow=rmax,ncol=cmax) Rnames=paste("AR",0:(ar.max),sep="-") Cnames=paste("MA",0:(ma.max),sep="-") rownames(corref)=Rnames colnames(corref)=Cnames rownames(cmj.table)=Rnames colnames(cmj.table)=Cnames rownames(pv)=Rnames colnames(pv)=Cnames rownames(mark)=Rnames colnames(mark)=Cnames for (m in 0:ar.max) { m1=m+1 for (j in 0:ma.max) { j1=j+1 if (m == 0 && j != 0) { racf=acf(z,plot=FALSE)$acf[1:(j+1)] lamb=racf[j+1]^2 corref[m1,j]=round(lamb,4) dmj=1 + 2*sum(racf[1:j]^2) cmj=-1*(n-m-j)*log(1.0 - lamb/dmj) pvalue =pchisq(cmj,1,lower.tail=FALSE) pv[m1,j]=round(pvalue,4) cmj.table[m1,j]=round(cmj,4) mark[m1,j]=ifelse(pvalue > alpha,"O","X") } else if (m != 0 && j == 0) { racf=pacf(z,plot=FALSE)$acf[1:(m+1)] lamb=racf[m+1]^2 corref[m1,j1]=round(lamb,4) dmj = 1 cmj=-1*(n-m-j)*log(1.0 - lamb/dmj) pvalue =pchisq(cmj,1,lower.tail=FALSE) pv[m1,j1]=round(pvalue,4) cmj.table[m1,j1]=round(cmj,4) mark[m1,j1]=ifelse(pvalue > alpha,"O","X") } else { mat1=matrix(0,nrow=m1,ncol=m1) mat2=matrix(0,nrow=m1,ncol=m1) mat3=matrix(0,nrow=m1,ncol=m1) mat4=matrix(0,nrow=m1,ncol=m1) for (t in (j+m+2):n) { tj1=t-j-1 ym1=ym(z,tj1,m) ym2=ym(z,t,m) mat1=mat1 + as.matrix(ym1)%*%ym1 mat2=mat2 + as.matrix(ym1)%*%ym2 mat3=mat3 + as.matrix(ym2)%*%ym2 mat4=mat4 + as.matrix(ym2)%*%ym1 } b1=solve(mat1)%*%mat2 b2=solve(mat3)%*%mat4 A=b2%*%b1 eig <-eigen(A) eig.val <-eig$values eig.val=Re(eig.val) eig.len=length(eig.val) eig.vector=eig$vectors lamb=min(eig.val) eig.vector0=eig.vector[,which.min(eig.val)] eig.vector0 = eig.vector0/eig.vector0[1] resid=(1:n)*0 for (t in (j+m+1):n) { z0=z[seq(t,t-m,-1)] resid[t]=sum(z0 * eig.vector0) } jm1=j + m + 1 rx=Re(resid[jm1:n]) racf=acf(rx,plot=FALSE)$acf[1:j] dmj=1 + 2*sum(racf^2) cmj=-1*(n-m-j)*log(1.0 - lamb/dmj) pvalue =pchisq(cmj,df=1,lower.tail=FALSE) corref[m1,j1]=round(lamb,4) pv[m1,j1]=round(pvalue,4) cmj.table[m1,j1]=round(cmj,4) mark[m1,j1]=ifelse(pvalue > alpha,"O","X") } } } cat("\n\nSCAN: Smallest CANonical Correlation Method for ARIMA(p,d,q)\n\n") cat("Estimates of Squared Canonical Correlation \n\n") print(corref) cat("\n\nC(m,j)\n\n") print(cmj.table) cat("\n\nChi-Square(1) Test p-value\n\n") print(pv) cat("\nSCAN Matrix \n\n") print(mark) plot(0:1,0:1,col="white",xlim=c(0,nrow(pv)-1),ylim=c(0,ncol(pv)-1),axes=FALSE,xlab="AR",ylab="MA") axis(1); axis(2) library(RColorBrewer) CL=brewer.pal(6, "RdBu")[c(1,2,3,5)] cpv=matrix(as.numeric(cut(as.vector(pv),c(-1,.01,.05,.1,2))),nrow(pv),ncol(pv)) for(i in 1:nrow(pv)){ for(j in 1:ncol(pv)){ polygon(c(i-1,i-1,i,i)-.5,c(j-1,j,j,j-1)-.5, col=CL[cpv[i,j]]) }} }

Consider the following simulated time series,

> s=arima.sim(n=200,model=list(ar=c(0,0,0,.4,0,0,0,.5),ma=c(0,0,1))) > plot(s,type="l")

The output is here

> arma.scan(s,6,6) SCAN: Smallest CANonical Correlation Method for ARIMA(p,d,q) Estimates of Squared Canonical Correlation MA-0 MA-1 MA-2 MA-3 MA-4 MA-5 MA-6 AR-0 0.0614 0.0104 0.1862 0.3516 0.0971 0.0128 0.0000 AR-1 0.0302 0.0294 0.1501 0.0943 0.0855 0.0127 0.0385 AR-2 0.3070 0.2781 0.2140 0.0006 0.1589 0.1884 0.2243 AR-3 0.1627 0.0037 0.1927 0.2311 0.1379 0.0207 0.0376 AR-4 0.2087 0.3947 0.3653 0.3075 0.1502 0.1364 0.1013 AR-5 0.1677 0.1219 0.0110 0.0263 0.0332 0.0350 0.0044 AR-6 0.0114 0.0485 0.0561 0.0427 0.0009 0.0089 0.0308 C(m,j) MA-0 MA-1 MA-2 MA-3 MA-4 MA-5 MA-6 AR-0 4.1161 0.6585 12.0315 20.6512 4.5388 0.5620 0.0000 AR-1 6.1127 1.9499 9.9356 4.9145 4.7219 0.4642 1.9015 AR-2 72.6011 19.1679 14.3512 0.0337 7.9668 9.6479 11.4573 AR-3 34.9724 0.2386 10.1620 13.4082 6.7875 0.8725 1.4071 AR-4 45.8691 27.5070 19.1422 20.2835 7.3339 5.5374 3.5874 AR-5 35.7981 8.0498 0.6280 1.3543 1.8470 1.7930 0.2338 AR-6 2.2147 3.1466 3.5990 1.9904 0.0511 0.4816 1.6440 Chi-Square(1) Test p-value MA-0 MA-1 MA-2 MA-3 MA-4 MA-5 MA-6 AR-0 0.0425 0.4171 0.0005 0.0000 0.0331 0.4534 0.0000 AR-1 0.0134 0.1626 0.0016 0.0266 0.0298 0.4957 0.1679 AR-2 0.0000 0.0000 0.0002 0.8543 0.0048 0.0019 0.0007 AR-3 0.0000 0.6252 0.0014 0.0003 0.0092 0.3503 0.2355 AR-4 0.0000 0.0000 0.0000 0.0000 0.0068 0.0186 0.0582 AR-5 0.0000 0.0046 0.4281 0.2445 0.1741 0.1806 0.6287 AR-6 0.1367 0.0761 0.0578 0.1583 0.8212 0.4877 0.1998 SCAN Matrix MA-0 MA-1 MA-2 MA-3 MA-4 MA-5 MA-6 AR-0 "O" "O" "X" "X" "O" "O" "X" AR-1 "O" "O" "X" "O" "O" "O" "O" AR-2 "X" "X" "X" "O" "X" "X" "X" AR-3 "X" "O" "X" "X" "X" "O" "O" AR-4 "X" "X" "X" "X" "X" "O" "O" AR-5 "X" "X" "O" "O" "O" "O" "O" AR-6 "O" "O" "O" "O" "O" "O" "O"

with the following graph

Of course, it is possible to ask for larger values,

> arma.scan(s,12,12)

The graph is now

# Temperatures Series as Random Walks

Last year, I did mention in a post that unit-root tests are dangerous, because they might lead us to strange models. For instance, in a post, I did obtain that the temperature observed in January 2013, in Montréal, might be considered as a random walk process (or at leat an integrated process). The code to extract the data has changed (since the website has been updated), so here, we use

library(RCurl) library(XML) options(RCurlOptions = list(useragent = "R")) HEURE=0:23 extracttemp=function(Y,M,D){ url=paste( "http://climate.weather.gc.ca/climateData/hourlydata_e.html?timeframe=1&Prov=QC&StationID=5415&Year=",Y,"&Month=", M,"&Day=",D,sep="") wp <- getURLContent(url) doc <- htmlParse(wp, asText = TRUE) docName(doc) <- url tmp <- readHTMLTable(doc) basejour=data.frame(Year=Y,Month=M,Day=D, Hour=HEURE,Temp=as.numeric(as.character(data.frame(tmp[2])[,2]))[2:25]) return(basejour)} B=NULL for(y in 1955:2013){ for(d in 1:31){ B=rbind(B,extracttemp(y,1,d))}}

Here are all the temperatures observed, and 2013,

plot(B$X,B$Temp,cex=.5,col="light blue",xlab="January, in Montreal",ylab="Temperature (Celsius)") I=which(B$Year==2013) lines(B$X[I],B$Temp[I],col="red")

In the previous post, one test only was used, and one year was considered. I was wondering if this behavior was observed only with temperature of 2013 (or not), and how the other tests (mentioned in a previous post too) were performing.

I might need a function, because those tests cannot be used if there is a missing value, even only one. So I did use the value observed one hour before (just to make sure that the tests can be done)

correcty=function(Y){ I=which(is.na(Y)) if(length(I)==0){Yc=Y} if(length(I)>0){Yc=Y;for(i in I) Yc[i]=Yc[i-1]} return(Yc) }

Now, we can compute the *p*-values, for all the years, and the three different three (keeping in mind that two test if the series is non-stationary, and one if the series is stationary)

DF=matrix(NA,2013-1954,3) library(urca) for(y in 1955:2013){ Z=B$Temp[which(B$Year==y)] Zc=correcty(Z) DF[y-1954,2]=as.numeric(pp.test(Zc)$p.value) DF[y-1954,1]=as.numeric(kpss.test(Zc)$p.value) DF[y-1954,3]=as.numeric(adf.test(Zc)$p.value) }

Visually, if red means stationary, and blue means non-stationary, we get

DFP=DF DFP[,1]=DF[,1]<.05 DFP[,2:3]=DF[,2:3]>.05 library(RColorBrewer) CL=brewer.pal(6, "RdBu") plot(0:1,0:1,xlim=c(1950,2015),ylim=c(0,3),axes=FALSE,xlab="",ylab="") axis(1) text(1952,.5,"KPSS") text(1952,1.5,"PP") text(1952,2.5,"ADF") for(y in 1955:2013){ for(i in 1:3){ polygon(y+c(-1,-1,1,1)/2.2,i-.5+c(-1,1,1,-1)/2.2,col=CL[1+(DFP[y-1954,i]==1)*5],border=NA)}}

Quite frequently, we conclude that the temperature is a random walk. Which does not make sense (from a physical point of view). But again, it might come from the fact that temperature are stationary, but with some fractional behavior (as suggested in the previous post).

# Unit Root Tests

This week, in the MAT8181 Time Series course, we’ve discussed unit root tests. According to Wold’s theorem, if is (weakly) stationnary then

where is the innovation process, and where is some deterministic series (just to get a result as general as possible). Observe that

as discussed in a previous post. To go one step further, there is also the Beveridge-Nelson decomposition : an integrated of order one process, defined as

can be represented as

a linear trend a random walk a stationary remaining term

i.e.

where is the polynomial with terms , where

For unit-root tests, we will use various representation of the process. In order to illustrate the implementation of those tests, consider the following series

> E=rnorm(240) > X=cumsum(E) > plot(X,type="l")

**Dickey Fuller (standard)**

Here, for the simple version of the Dickey-Fuller test, we assume that

and we would like to test if (or not). We can write the previous representation as

so we simply have to test if the regression coefficient in the linear regression is – or not – null. Which can be done with Student’s test. If we consider the previous model without the linear drift, we have to consider the following regression

> lags=0 > z=diff(X) > n=length(z) > z.diff=embed(z, lags+1)[,1] > z.lag.1=X[(lags+1):n] > summary(lm(z.diff~0+z.lag.1 )) Call: lm(formula = z.diff ~ 0 + z.lag.1) Residuals: Min 1Q Median 3Q Max -2.84466 -0.55723 -0.00494 0.63816 2.54352 Coefficients: Estimate Std. Error t value Pr(>|t|) z.lag.1 -0.005609 0.007319 -0.766 0.444 Residual standard error: 0.963 on 238 degrees of freedom Multiple R-squared: 0.002461, Adjusted R-squared: -0.00173 F-statistic: 0.5873 on 1 and 238 DF, p-value: 0.4442

Our testing procedure will be based on the Student’s *t* value,

> summary(lm(z.diff~0+z.lag.1 ))$coefficients[1,3] [1] -0.7663308

which is exactly the value computed using

> library(urca) > df=ur.df(X,type="none",lags=0) > df ############################################################### # Augmented Dickey-Fuller Test Unit Root / Cointegration Test # ############################################################### The value of the test statistic is: -0.7663

The interpretation of this value can be done using critical values (99%, 95%, 90%)

> qnorm(c(.01,.05,.1)/2) [1] -2.575829 -1.959964 -1.644854

If the statistics exceeds those values, then the series is not stationnary, since we cannot reject the assumption that . So we might conclude that there is a unit root. Actually, those critical values are obtained using

> summary(df) ############################################### # Augmented Dickey-Fuller Test Unit Root Test # ############################################### Test regression none Call: lm(formula = z.diff ~ z.lag.1 - 1) Residuals: Min 1Q Median 3Q Max -2.84466 -0.55723 -0.00494 0.63816 2.54352 Coefficients: Estimate Std. Error t value Pr(>|t|) z.lag.1 -0.005609 0.007319 -0.766 0.444 Residual standard error: 0.963 on 238 degrees of freedom Multiple R-squared: 0.002461, Adjusted R-squared: -0.00173 F-statistic: 0.5873 on 1 and 238 DF, p-value: 0.4442 Value of test-statistic is: -0.7663 Critical values for test statistics: 1pct 5pct 10pct tau1 -2.58 -1.95 -1.62

The problem with R is that there are several packages that can be used for unit root tests. Just to mention another one,

> library(tseries) > adf.test(X,k=0) Augmented Dickey-Fuller Test data: X Dickey-Fuller = -2.0433, Lag order = 0, p-value = 0.5576 alternative hypothesis: stationary

We do have here also a test where the null hypothesis is that there is a unit root. But the *p*-value is quite different. What is odd is that we have

> 1-adf.test(X,k=0)$p.value [1] 0.4423705 > df@testreg$coefficients[4] [1] 0.4442389

(but I think it is a coincidence).

**Augmented Dickey Fuller**

It is possible to had some lags in the regression. For instance, we can consider

Again, we have to check if one coefficient is null, or not. And this can be done using Student’s* t* test.

> lags=1 > z=diff(X) > n=length(z) > z.diff=embed(z, lags+1)[,1] > z.lag.1=X[(lags+1):n] > k=lags+1 > z.diff.lag = embed(z, lags+1)[, 2:k] > summary(lm(z.diff~0+z.lag.1+z.diff.lag )) Call: lm(formula = z.diff ~ 0 + z.lag.1 + z.diff.lag) Residuals: Min 1Q Median 3Q Max -2.87492 -0.53977 -0.00688 0.64481 2.47556 Coefficients: Estimate Std. Error t value Pr(>|t|) z.lag.1 -0.005394 0.007361 -0.733 0.464 z.diff.lag -0.028972 0.065113 -0.445 0.657 Residual standard error: 0.9666 on 236 degrees of freedom Multiple R-squared: 0.003292, Adjusted R-squared: -0.005155 F-statistic: 0.3898 on 2 and 236 DF, p-value: 0.6777 > summary(lm(z.diff~0+z.lag.1+z.diff.lag ))$coefficients[1,3] [1] -0.7328138

This value is the one obtained using

> df=ur.df(X,type="none",lags=1) > summary(df) ############################################### # Augmented Dickey-Fuller Test Unit Root Test # ############################################### Test regression none Call: lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag) Residuals: Min 1Q Median 3Q Max -2.87492 -0.53977 -0.00688 0.64481 2.47556 Coefficients: Estimate Std. Error t value Pr(>|t|) z.lag.1 -0.005394 0.007361 -0.733 0.464 z.diff.lag -0.028972 0.065113 -0.445 0.657 Residual standard error: 0.9666 on 236 degrees of freedom Multiple R-squared: 0.003292, Adjusted R-squared: -0.005155 F-statistic: 0.3898 on 2 and 236 DF, p-value: 0.6777 Value of test-statistic is: -0.7328 Critical values for test statistics: 1pct 5pct 10pct tau1 -2.58 -1.95 -1.62

And again, other pckages can be used:

> adf.test(X,k=1) Augmented Dickey-Fuller Test data: X Dickey-Fuller = -1.9828, Lag order = 1, p-value = 0.5831 alternative hypothesis: stationary

Hopefully, the conclusion is the same (we should reject the assumption that the series is stationary, but I am not sure about the computation of the *p*-value).

**Augmented Dickey Fuller with trend and drift**

So far, we have not included the drift in our model. But this is simple to do (this will be called the augmented version of the previous procedure): we just have to include a constant in the regression,

> summary(lm(z.diff~1+z.lag.1+z.diff.lag )) Call: lm(formula = z.diff ~ 1 + z.lag.1 + z.diff.lag) Residuals: Min 1Q Median 3Q Max -2.91930 -0.56731 -0.00548 0.62932 2.45178 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.29175 0.13153 2.218 0.0275 * z.lag.1 -0.03559 0.01545 -2.304 0.0221 * z.diff.lag -0.01976 0.06471 -0.305 0.7603 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.9586 on 235 degrees of freedom Multiple R-squared: 0.02313, Adjusted R-squared: 0.01482 F-statistic: 2.782 on 2 and 235 DF, p-value: 0.06393

The statistics of interest are obtained here considering some analysis of variance outputs, where this model is compared with the one without the integrated part, and the drift,

> summary(lm(z.diff~1+z.lag.1+z.diff.lag ))$coefficients[2,3] [1] -2.303948 > anova(lm(z.diff ~ z.lag.1 + 1 + z.diff.lag),lm(z.diff ~ 0 + z.diff.lag))$F[2] [1] 2.732912

Those two values are the ones obtained also with

> df=ur.df(X,type="drift",lags=1) > summary(df) ############################################### # Augmented Dickey-Fuller Test Unit Root Test # ############################################### Test regression drift Call: lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag) Residuals: Min 1Q Median 3Q Max -2.91930 -0.56731 -0.00548 0.62932 2.45178 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.29175 0.13153 2.218 0.0275 * z.lag.1 -0.03559 0.01545 -2.304 0.0221 * z.diff.lag -0.01976 0.06471 -0.305 0.7603 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.9586 on 235 degrees of freedom Multiple R-squared: 0.02313, Adjusted R-squared: 0.01482 F-statistic: 2.782 on 2 and 235 DF, p-value: 0.06393 Value of test-statistic is: -2.3039 2.7329 Critical values for test statistics: 1pct 5pct 10pct tau2 -3.46 -2.88 -2.57 phi1 6.52 4.63 3.81

And we can also include a linear trend,

> temps=(lags+1):n > summary(lm(z.diff~1+temps+z.lag.1+z.diff.lag )) Call: lm(formula = z.diff ~ 1 + temps + z.lag.1 + z.diff.lag) Residuals: Min 1Q Median 3Q Max -2.87727 -0.58802 -0.00175 0.60359 2.47789 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.3227245 0.1502083 2.149 0.0327 * temps -0.0004194 0.0009767 -0.429 0.6680 z.lag.1 -0.0329780 0.0166319 -1.983 0.0486 * z.diff.lag -0.0230547 0.0652767 -0.353 0.7243 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.9603 on 234 degrees of freedom Multiple R-squared: 0.0239, Adjusted R-squared: 0.01139 F-statistic: 1.91 on 3 and 234 DF, p-value: 0.1287 > summary(lm(z.diff~1+temps+z.lag.1+z.diff.lag ))$coefficients[3,3] [1] -1.98282 > anova(lm(z.diff ~ z.lag.1 + 1 + temps+ z.diff.lag),lm(z.diff ~ 1+ z.diff.lag))$F[2] [1] 2.737086

while R function returns

> df=ur.df(X,type="trend",lags=1) > summary(df) ############################################### # Augmented Dickey-Fuller Test Unit Root Test # ############################################### Test regression trend Call: lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag) Residuals: Min 1Q Median 3Q Max -2.87727 -0.58802 -0.00175 0.60359 2.47789 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.3227245 0.1502083 2.149 0.0327 * z.lag.1 -0.0329780 0.0166319 -1.983 0.0486 * tt -0.0004194 0.0009767 -0.429 0.6680 z.diff.lag -0.0230547 0.0652767 -0.353 0.7243 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.9603 on 234 degrees of freedom Multiple R-squared: 0.0239, Adjusted R-squared: 0.01139 F-statistic: 1.91 on 3 and 234 DF, p-value: 0.1287 Value of test-statistic is: -1.9828 1.8771 2.7371 Critical values for test statistics: 1pct 5pct 10pct tau3 -3.99 -3.43 -3.13 phi2 6.22 4.75 4.07 phi3 8.43 6.49 5.47

**KPSS test**

Here, in the KPSS testing procedure, two models can be considerd : with a drift, or with a linear trend. Here, the null hypothesis is that the series is stationnary.

With a drift, the code is

> summary(ur.kpss(X,type="mu")) ####################### # KPSS Unit Root Test # ####################### Test is of type: mu with 4 lags. Value of test-statistic is: 0.972 Critical value for a significance level of: 10pct 5pct 2.5pct 1pct critical values 0.347 0.463 0.574 0.73

while it will be, in the case there is a trend

> summary(ur.kpss(X,type="tau")) ####################### # KPSS Unit Root Test # ####################### Test is of type: tau with 4 lags. Value of test-statistic is: 0.5057 Critical value for a significance level of: 10pct 5pct 2.5pct 1pct critical values 0.119 0.146 0.176 0.216

One more time, it is possible to use another package to get the same test (but again, a different output)

> kpss.test(X,"Level") KPSS Test for Level Stationarity data: X KPSS Level = 1.1997, Truncation lag parameter = 3, p-value = 0.01 > kpss.test(X,"Trend") KPSS Test for Trend Stationarity data: X KPSS Trend = 0.6234, Truncation lag parameter = 3, p-value = 0.01

At least, there is some kind of consistency, since we keep rejecting the stationnary assumption, for that series.

**Philipps-Perron test**

The Philipps-Perron test is based on the ADF procedure. The code is here

> PP.test(X) Phillips-Perron Unit Root Test data: X Dickey-Fuller = -2.0116, Truncation lag parameter = 4, p-value = 0.571

with again, a possible alternative with the other package

> pp.test(X) Phillips-Perron Unit Root Test data: X Dickey-Fuller Z(alpha) = -7.7345, Truncation lag parameter = 4, p-value = 0.6757 alternative hypothesis: stationary

**Comparison**

I will not spend more time comparing the different codes, in R, to run those tests. Let us spend some additional time on a quick comparison of those three procedure. Let us generate some autoregressive processes, with more or less autocorrelation, as well as some random walk, and let us see how those tests perform :

> n=100 > AR=seq(1,.7,by=-.01) > P=matrix(NA,3,31) > M1=matrix(NA,1000,length(AR)) > M2=matrix(NA,1000,length(AR)) > M3=matrix(NA,1000,length(AR)) > for(i in 1:(length(AR)+1)){ + for(s in 1:1000){ + if(i==1) X=cumsum(rnorm(n)) + if(i!=1) X=arima.sim(n=n,list(ar=AR[i])) + library(urca) + M2[s,i]=as.numeric(pp.test(X)$p.value) + M1[s,i]=as.numeric(kpss.test(X)$p.value) + M3[s,i]=as.numeric(adf.test(X)$p.value) + }}

Here, we would like to count how many times the *p*-value of our tests exceed 5%,

> prop05=function(x) mean(x>.05) + P[1,]=1-apply(M1,2,prop05) + P[2,]=apply(M2,2,prop05) + P[3,]=apply(M3,2,prop05) + } > plot(AR,P[1,],type="l",col="red",ylim=c(0,1),ylab="proportion of non-stationnary + series",xlab="autocorrelation coefficient") > lines(AR,P[2,],type="l",col="blue") > lines(AR,P[3,],type="l",col="green") > legend(.7,1,c("ADF","KPSS","PP"),col=c("green","red","blue"),lty=1,lwd=1)

We can see here how poorly Dickey-Fuller test behave, since a 50% (at least) of our autoregressive processes are considered as non-stationnary.

# Inference for ARMA(p,q) Time Series

As we mentioned in our previous post, as soon as we have a moving average part, inference becomes more complicated. Again, to illustrate, we do not need a two general model. Consider, here, some process,

where is some white noise, and assume further that .

> theta=.7 > phi=.5 > n=1000 > Z=rep(0,n) > set.seed(1) > e=rnorm(n) > for(t in 2:n) Z[t]=phi*Z[t-1]+e[t]+theta*e[t-1] > Z=Z[800:1000] > plot(Z,type="l")

- A
**two step procedure**

To start with something simple, assume that we did miss the moving average component,

The estimator of – by least squares – is not longer consistent. But still. We can still compute it

> base=data.frame(Y=Z[2:n],X=Z[1:(n-1)]) > regression=lm(Y~0+X,data=base) > summary(regression) Call: lm(formula = Y ~ 0 + X, data = base) Residuals: Min 1Q Median 3Q Max -3.2445 -0.7909 0.0626 0.9707 3.0685 Coefficients: Estimate Std. Error t value Pr(>|t|) X 0.69571 0.05101 13.64 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.225 on 199 degrees of freedom (799 observations deleted due to missingness) Multiple R-squared: 0.4832, Adjusted R-squared: 0.4806 F-statistic: 186 on 1 and 199 DF, p-value: < 2.2e-16

and then, we cancompute the autocorrelation of the noise,

> n=200 > cor(residuals(regression)[2:n],residuals(regression)[1:(n-1)]) [1] 0.2663076

or more formally, use Durbin-Watson estimator, to get autocorrelation of the noise (and some significance test)

> library(car) > durbinWatsonTest(regression) lag Autocorrelation D-W Statistic p-value 1 0.2656555 1.46323 0 Alternative hypothesis: rho != 0

The point, here, is that we would like to assume that

meaning that should be some process. And

i.e. is a root of this quadratic problem,

> polyroot(c(1,-1/cor(residuals(regression)[2:n],residuals(regression)[1:(n-1)]),1)) [1] 0.2884681+0i 3.4665883+0i

Here, we do have two positive roots. I would go for the one smaller than one, in order to be able to invert the polynomial, if necessary…

- Use of the
**empirical autocorrelation function**

An alternative might be to use properties of the autocorrelation function,

and

Again, we have a set of three equations, with three unknown parameters. Numerically, it is possible to find some roots. If we run the code, we get

> v=c(as.numeric(acf(Z)$acf[2:3]),1)*var(Z) > library(rootSolve) > seteq=function(x){ + F1=v[1]-x[3]^2*(x[2]^2+2*x[1]*x[2]+1)/(1-x[1]^2) + F2=v[2]-(x[1]*v[1]+x[2]*x[3]^2) + F3=v[3]-x[1]*v[2] + return(c(F1,F2,F3))} > multiroot(f=seteq,start=c(.1,.1,1)) $root [1] 3.643734 -3.188145 1.427759 $f.root [1] 1.371170e-11 -3.714573e-11 0.000000e+00 $iter [1] 8 $estim.precis [1] 1.695248e-11

Here, we have a situation…

- Use of
**least square techniques**

We can use, here, the algorithm described in the context of processes.

> V=function(p){ + phi=p[1] + theta=p[2] + u=rep(0,length(Z)) + for(t in 2:length(Z)) u[t]=Z[t]-phi*Z[t-1]-theta*u[t-1] + return(sum(u^2)) + } > p=optim(par=c(.1,.1),V)$par [1] 0.3637783 0.7773845 > coef=c(p,sqrt(V(p)/(length(Z))))

which is not so bad. Actually, if we run that procecure on 1,000 samples, we get the following output

- Use of
**maximum likelihood techniques**

Last, but not least, one more time, we can use (global) maximum likelihood techniques, since the process is a Gaussian process (all finite dimensional vector will have a joint Gaussian distribution) if we assume that the noise is Gaussian.

> library(mnormt) > GlobalLogLik=function(A,TS){ + n=length(TS) + phi=A[1]; theta=A[2] + sigma=A[3] + SIG=matrix(0,n,n) + rho=rep(0,n) + rho[1]=sigma^2*(theta^2+2*phi*theta+1)/(1-phi^2) + rho[2]=phi*rho[1]+theta*sigma^2 + for(h in 3:n) rho[h]=phi*rho[h-1] + for(i in 1:n){for(j in 1:n){ + SIG[i,j]=rho[abs(i-j)+1]}} + return(dmnorm(TS,rep(0,n),SIG,log=TRUE))} > LogL=function(A) -GlobalLogLik(A,TS=Z) > optim(c(.1,.1,1),LogL)$par [1] 0.3890991 0.7672036 1.0731340

It works fine, one more time. But maybe we got lucky here. We’ve seen in the post on autoregressive time series that the algorithm might fell if the time series is not stationary. In order to avoid such problems, we can consider a constraint optimization problem, where we simply recall that ,

> U=matrix(c(1,-1,0,0,0,0),2,3) > C=-c(.999,.999) > constrOptim(c(.1,.1,1),LogL,grad=NULL,ui=U,ci=C) $par [1] 0.3890991 0.7672036 1.0731340 $value [1] 300.1956 $counts function gradient 118 NA $convergence [1] 0 $message NULL $outer.iterations [1] 2 $barrier.value [1] -1.536358e-05

If we run that algorithm 1,000 times, on simulated time series (with the same parameters), we get

# Inference for MA(q) Time Series

Yesterday, we’ve seen how inference for time series was possible. I started with that one because it is actually the simple case. For instance, we can use ordinary least squares. There might be some possible bias (see e.g. White (1961)), but asymptotically, estimators are fine (consistent, with asymptotic normality). But when the noise is (auto)correlated, then it is more complex. So, consider here some time series

for some white noise .

> theta1=.25 > theta2=.7 > n=1000 > set.seed(1) > e=rnorm(n) > Z=rep(0,n) > for(t in 3:n) Z[t]=e[t]+theta1*e[t-1]+theta2*e[t-2] > Z=Z[800:1000] > plot(Z,type="l")

- Using the
**empirical autocorrelations**

The first idea might be to use the first two (empirical) autocorrelations (the two that are supposed to be – theoretically – non null).

with when . We also have the following relationship on the variance of the process

With those three equations, for three unknown parameters, , and , we *simply* have to solve (numerically) that system of equations,

> v=c(as.numeric(acf(Z)$acf[2:3]),var(Z)) > v [1] 0.1658760 0.3823053 1.6379498 > library(rootSolve) > seteq=function(x){ + F1=v[1]-(x[1]+x[1]*x[2])/(1+x[1]^2+x[2]^2) + F2=v[2]-(x[2])/(1+x[1]^2+x[2]^2) + F3=v[3]-(1+x[1]^2+x[2]^2)*x[3]^2 + return(c(F1,F2,F3))} > multiroot(f=seteq,start=c(.1,.1,1)) $root [1] 0.1400579 0.4766699 1.1461636 $f.root [1] 7.876355e-10 4.188458e-09 -2.839977e-09 $iter [1] 5 $estim.precis [1] 2.605357e-09

We are a bit far away from the *true* values, used to generate our sample. And if we consider 1,000 sample (instead of only one), we still have the bias, and a large variance for our three estimators,

- Using
**least square techniques**

We can try something quite different here. The problem we have is that we do not observe the noise , we only observe our series . But we can try to rebuild that series (call it since we’re not sure it will be a reconstruction of the noise). As suggested in Box & Jenkins (1967), assume that the first two values are null. And then, use

and then, we can use least square techniques

The code will be

> V=function(p){ + theta1=p[1] + theta2=p[2] + u=rep(0,length(Z)) + for(t in 3:length(Z)) u[t]=Z[t]-theta1*u[t-1]-theta2*u[t-2] + return(sum(u^2)) + }

If we try to minimize the sum of the squares of the residuals, we get

> optim(par=c(.1,.1),V) $par [1] 0.2751667 0.6723909 $value [1] 225.8104 $counts function gradient 77 NA $convergence [1] 0 $message NULL

which is close to the *true* value. Another good thing is that, if we compare that rebuilt noise with the true one (since we actually have it), then we have the same vector,

> plot(e[800:1000],col="blue",type="l") > theta1=0.2751667 > theta2=0.6723909 > u=rep(0,length(Z)) > for(t in 3:length(Z)) u[t]=Z[t]-theta1*u[t-1]-theta2*u[t-2] > lines(1:201,u,col="red")

So far, so good. And if we look at 1,000 samples, we get

It looks like we have some bias here. And since the two estimators should be negatively correlated, one over-estimates, while the other one under-estimates.

- Using the
**(global) maximum likelihood technique**

And a final method might be to use the maximum likelihood technique (globally). Again, if we assume that we have a Gaussian i.i.d noise, then the vector is Gaussian, with a simple variance matrix (since a lot of elements will be null),

> library(mnormt) > GlobalLogLik=function(A,TS){ + n=length(TS) + theta1=A[1]; theta2=A[2] + sigma=A[3] + SIG=matrix(0,n,n) + rho=rep(0,n) + rho[1]=1 + rho[2]=(theta1+theta1*theta2)/(1+theta1^2+theta2^2) + rho[3]=(theta2)/(1+theta1^2+theta2^2) + for(i in 1:n){for(j in 1:n){ + SIG[i,j]=rho[abs(i-j)+1]}} + gamma0=(1+theta1^2+theta2^2)*sigma^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) $par [1] 0.2584144 0.6826530 1.0669820 $value [1] 298.8699 $counts function gradient 86 NA $convergence [1] 0 $message NULL

Here, the values that minimize the likelihood are rather close to the ones used to generate our sample. And if we run this algorithm on 1,000 samples, we can see that those estimates are fine,

I could not find other ideas, to estimate those parameters. I guess we can use the partial autocorrelation function, since we have relationships that can be related to Yule-Walker equations for time series.

# 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] 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[1]*base$X1+PHI[2]*base$X2) > sd(estWN) [1] 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[1]; phi2=A[2] + sigma=A[3]; 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 [1] 0.2339589 0.6285002 1.0565613 $value [1] 293.3042 $counts function gradient 106 NA $convergence [1] 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[1]; phi2=A[2] + sigma=A[3] + SIG=matrix(0,n,n) + rho=rep(0,n) + rho[1]=1 + rho[2]=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) > constrOptim(c(.1,.1,1),LogL,grad=NULL,ui=U,ci=C) $par [1] 0.2238892 0.6342850 1.0613388 $value [1] 297.9202 $counts function gradient 108 NA $convergence [1] 0 $message NULL $outer.iterations [1] 2 $barrier.value [1] 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).

# Causal Autoregressive Time Series

In the MAT8181 graduate course on Time Series, we will discuss (almost) only **causal** models. For instance, with ,

with some white noise , those models are obtained when . In that case, we’ve seen that was actually the innovation process, and we can write

which is actually a mean-square convergent series (using simple Analysis arguments on series). From that expression, we can easily see that is stationary, since (which does not depend on ) and

(which does not depend on ).

Consider now the case where . Clearly, we have some problem here, since

cannot be defined (the series does not converge, in ). Nevertheless, it is still possible to write

But it is possible to iterate (as in the previous case) and write

which is actually well defined. And in that case, the sequence of random variables obtained from this equation is the unique stationary solution of the recursive equation . This might be confusing, but the thing is this solution should not be confused with the *usual* non-stationary solution of obtained from . As in the code writen to generate a time series, *from some starting value *in the previous post.

Now, let us spent some time with this stationary time series, considered as *unatural *in Brockwell and Davis (1991). One point is that, in the previous case (where ) was the innovation process. So variable was not correlated with the future of the noise, . Which is not the case when .

All that looks nice, if you’re willing to understand thing at some theoretical level. What does all that mean from a computational perspective ? Consider some white noise (this noise actually does exist whatever you want to define, based on that time series)

> n=10000 > e=rnorm(n) > plot(e,type="l",col="red")

If we look at the simple case, to start with,

> phi=.8 > X=rep(0,n) > for(t in 2:n) X[t]=phi*X[t-1]+e[t]

The time series – the latest 1,000 observations – looks like

Now, if we use the cumulated sum of the noise,

> Y=rep(0,n) > for(t in 2:n) Y[t]=sum(phi^((0:(t-1)))*e[t-(0:(t-1))])

we get

Which is exactly the same process ! This should not surprise us because that’s what the theory told us. Now, consider the problematic case, where

> phi=1.1 > X=rep(0,n) > for(t in 2:n) X[t]=phi*X[t-1]+e[t]

Clearly, that series is non-stationary (just look at the first 1,000 values)

Now, if we look at the series obtained from the cumulated sum of *future* values of the noise

> Y=rep(0,n) > for(t in 1:(n-1)) Y[t]=sum((1/phi)^((1:(n-t)))*e[t+(1:(n-t))])

We get something which is, actually, stationary,

So, what is this series exactly ? If you look that the autocorrelation function,

> acf(Y)

we get the autocorrelation function of a (stationary) process,

> acf(Y)[1] Autocorrelations of series ‘Y’, by lag 1 0.908 > 1/phi [1] 0.9090909

Observe that there is a white noise – call it – such that

This is what we call the *canonical *form of the stationary process .

# 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)