Tag Archives: AR(1)

Forecasts with ARIMA Models

In our time series class this morning, I was discussing forecasts with ARIMA Models. Consider some simple stationnary AR(1) simulated time series

> n=95
> set.seed(1)
> E=rnorm(n)
> X=rep(0,n)
> phi=.85
> for(t in 2:n) X[t]=phi*X[t-1]+E[t]
> plot(X,type="l")

If we fit an AR(1) model,

> model=arima(X,order=c(1,0,0),
+             include.mean = FALSE)
> P=predict(model,n.ahead=20)
> plot(P$pred)
> lines(P$pred+2*P$se,col="red")
> lines(P$pred-2*P$se,col="red")
> abline(h=0,lty=2)
> abline(h=2*P$se[20],lty=2,col="red")
> abline(h=-2*P$se[20],lty=2,col="red")

we observe the exponential decay of the forecast towards 0, and the increasing confidence interval (where the variance increases, from the variance of the white noise to the variance of the stationnary time series). Plain lines are conditional forecast (given our latest observation, since the AR(1) is a first order Markov process), and dotted lines are unconditional. Let us store some values, to use them as benchmark

> s=P$se[20]
> y=P$pred

If we fit a MA(1) model

> model=arima(X,order=c(0,0,1),
+             include.mean = FALSE)
> P=predict(model,n.ahead=20)
> plot(P$pred)
> lines(P$pred+2*P$se,col="red")
> lines(P$pred-2*P$se,col="red")
> abline(h=0,lty=2)
> abline(h=2*s,lty=2,col="red")
> abline(h=-2*s,lty=2,col="red")
> lines(y,col="grey")

after two lags, the forecast is null, and the (conditional) variance remains constant. But if we consider a moving average process with a longer order,

> model=arima(X,order=c(0,0,14),
+             include.mean = FALSE)
> P=predict(model,n.ahead=20)
> plot(P$pred)
> lines(P$pred+2*P$se,col="red")
> lines(P$pred-2*P$se,col="red")
> abline(h=0,lty=2)
> abline(h=2*s,lty=2,col="red")
> abline(h=-2*s,lty=2,col="red")
> lines(y,col="grey")

we get an output that can be compared with the AR(1) processes. Which makes sense since our AR(1) process can also be seen as a MA(∞), with infinite order.

But if we think that our time series is not stationary, an we fit an integrated model

> model=arima(X,order=c(0,1,0),
+             include.mean = FALSE)
> P=predict(model,n.ahead=20)
> plot(P$pred)
> lines(P$pred+2*P$se,col="red")
> lines(P$pred-2*P$se,col="red")
> abline(h=0,lty=2)
> abline(h=2*s,lty=2,col="red")
> abline(h=-2*s,lty=2,col="red")
> lines(y,col="grey")

we observe the (standard) martingale property: the forecast is flat, and the confidence interval keeps increasing, and actually, the variance increases towards infinity (at a linear rate). So one should be very careful when differentiation a time series… it will have a huge impact on the forecasts….

Variance of the Average of a Sequence

In the case where  are i.i.d. random variables, then

Now, what if  are identically distributed, but no longer independent. What if we have an autoregressive process? Assume that


can be written

Here, we will express the variance as a function of  and , but it is possible to use also , since, in the context of an ,

Now, since  we get

which can be simplified, since


So, the variance of the mean can be writen as

Observe that if  is large enough,

This asymptotic relationship is well known actually. A simple way to get it is the following. One can can write

or equivalently

But actually, the first relationship is probably more interesting to get an asymptotic approximation,

In the context of an  process, this can be writen

Thus, we get the following well-known relationship

In the case where  is an i.i.d. sequence, i.e. , then we get the relationship mentioned initially. And in the case of a random walk… unfortunately, we cannot use that relationship. But observe that


which can be written

If we compare the true value and the approximation, we get the following graph,

> V=function(phi,s2=1,n=100){
+ g0=s2/(1-phi^2)
+ if(phi<1){
+ if(phi==0){v1=g0/n}
+ if(phi>0){v1=g0/n^2*(n+2*((n-1)*
+ phi^(-1)-n+phi^(n-1))/(phi^(-1)-1)^2)}
+ v2=g0/n*(1+phi)/(1-phi)
+ }
+ if(phi==1){
+ v1=(2*n+1)*(n+1)*s2/(6*n)
+ v2=NA
+ }
+ return(c(v1,v2))}
> Vphi=function(phi) V(phi,1,100)
> x=seq(.01,1,by=.02)
> M=matrix(unlist(lapply(x,V)),nrow=2)
> plot(x,M[1,],type="l",col="red",log="y",
+ ylab="Variance of the average (log scale)",
+ xlab="Autoregressive coefficient")
> lines(x,M[2,],col="blue")

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,


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


So finally, what we got here is


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.

Causal Autoregressive Time Series

In the MAT8181 graduate course on Time Series, we will discuss (almost) only causal models. For instance, with https://latex.codecogs.com/gif.latex?AR(1),


with some white noise https://latex.codecogs.com/gif.latex?(\varepsilon_t), those models are obtained when https://latex.codecogs.com/gif.latex?\vert%20\phi\vert%20%3C1. In that case, we’ve seen that https://latex.codecogs.com/gif.latex?(\varepsilon_t) 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 https://latex.codecogs.com/gif.latex?(X_t) is stationary, since https://latex.codecogs.com/gif.latex?\mathbb{E}(X_t)=0 (which does not depend on https://latex.codecogs.com/gif.latex?t) and

https://latex.codecogs.com/gif.latex?\text{cov}(X_t,X_{t-h})=\frac{\phi^h}{1-\phi^2}\sigma^2(which does not depend on https://latex.codecogs.com/gif.latex?t).

Consider now the case where https://latex.codecogs.com/gif.latex?\vert%20\phi\vert%20%3E1. Clearly, we have some problem here, since


cannot be defined (the series does not converge, in https://latex.codecogs.com/gif.latex?L^2). Nevertheless, it is still possible to write

https://latex.codecogs.com/gif.latex?X_t=\frac{1}{\phi}%20X_{t{\color{Red}%20+1}}{\color{Red}%20-\frac{1}{\phi}}\varepsilon_{t{\color{Red}%20+1}}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 https://latex.codecogs.com/gif.latex?(X_t) obtained from this equation is the unique stationary solution of the recursive equation https://latex.codecogs.com/gif.latex?X_t=\phi%20X_{t-1}+\varepsilon_t. This might be confusing, but the thing is this solution should not be confused with the usual non-stationary solution of https://latex.codecogs.com/gif.latex?X_t=\phi%20X_{t-1}+\varepsilon_t obtained from https://latex.codecogs.com/gif.latex?X_0. As in the code writen to generate a time series, from some starting value https://latex.codecogs.com/gif.latex?X_0 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 https://latex.codecogs.com/gif.latex?\vert%20\phi\vert%20%3C1) https://latex.codecogs.com/gif.latex?(\varepsilon_t) was the innovation process. So variable https://latex.codecogs.com/gif.latex?X_t was not correlated with the future of the noise, https://latex.codecogs.com/gif.latex?\sigma\{\varepsilon_{t+1},\varepsilon_{t+2},\cdots\}. Which is not the case when https://latex.codecogs.com/gif.latex?\vert%20\phi\vert%20%3E1.

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 https://latex.codecogs.com/gif.latex?\vert%20\phi\vert%20%3E1

> 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) https://latex.codecogs.com/gif.latex?AR(1) process,

> acf(Y)[1]

Autocorrelations of series ‘Y’, by lag


> 1/phi
[1] 0.9090909

Observe that there is a white noise – call it https://latex.codecogs.com/gif.latex?(\eta_t) – such that


This is what we call the canonical form of the stationary process https://latex.codecogs.com/gif.latex?(X_t).

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 https://latex.codecogs.com/gif.latex?AR(1) – 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))}


> graphar1(.8)


> graphar1(-.7)

(with, on the bottom, the root of the characteristic polynomial, the value of the parameter https://latex.codecogs.com/gif.latex?\phi_{1}, the autocorrelation function https://latex.codecogs.com/gif.latex?h\mapsto\rho(h) and the partial autocorrelation function https://latex.codecogs.com/gif.latex?h\mapsto\psi(h)).

Of course, it is possible to do something similar with https://latex.codecogs.com/gif.latex?AR(2) 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)


> graphar2(-1.4,-.7)

Unit root, or not ? is it a big deal ?

Consider a time series, generated using

for(t in 2:240){X[t]=rho*X[t-1]+E[t]}

The idea is to assume that an autoregressive model can be considered, but we don’t know the value of the parameter. More precisely, we can’t choose if the parameter is either one (and the series is integrated), or some value strictly smaller than 1 (and the series is stationary). Based on past observations, the higher the autocorrelation, the lower the variance of the noise.

rhoest=0.9; H=260

Now that we have a model, consider the following forecast, including a confidence interval,

Here, forecasts can be derived, with any kind of possible autoregressive coefficient, from 0.7 to 1. I.e. we can chose to model the time series either with a stationary, or an integrated series,

As we can see above, assuming either that the series is stationary (parameter lower – strictly – than 1) or integrated (parameter equal to 1), the shape of the prediction can be quite different. So yes, assuming an integrated model is a big deal, since it has a strong impact on predictions.

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

phi1= .85
for(t in 2:n){X[t]=phi1*X[t-1]+bruit[t]}

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
for(t in 3:n){

Histoire de changer un peu, on peut regarder un processus moyenne mobile au premier ordre,  où  est un paramètre dans .

theta1=  .8
for(t in 2:n){

ou une moyenne mobile du second ordre,

theta1= -.6
theta2=  .5
for(t in 3:n){


Trop de probabilités tue les probabilités

Hier matin, au moment de quitter la maison, j’ai jeté un œil à la météo (car à Montréal le temps est très changeant, pas seulement d’un jour sur l’autre, comme je l’avis mentionné en arrivant en septembre, ici). En voyant les gros nuages gris dans le ciel, en me réveillant, je me suis demandé s’il fallait mettre un manteau et des bottes en caoutchouc aux enfants (pour aller l’école). Sur le site de météo, on a en effet la probabilité de pluie, i.e. de P.D.P.

Si on suppose que les probabilités sont indépendantes d’heure en heure, la probabilité qu’il va pleuvoir dans la matinée (dans les 5 prochaines heures pour faire simple, ou plus concret), va s’écrire comme le complémentaire de la probabilité qu’il va faire beau toute la matinée, autrement dit, avec des probabilités de pluie de l’ordre de 25%, par heure, la probabilité qu’il pleuve pendant les prochaines heures est de 75% environ

> 1-((1-.25))^5
[1] 0.7626953

Bref, on est loin de 25%. Essayons d’aller un peu plus loin pour mieux comprendre…
Notons https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/correlation-pluie03.gif la variable indicatrice indiquant s’il pleut, ou pas (1 s’il pleut, et 0 sinon) à l’instant https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/correlation-pluie04.gif. Supposons qu’il existe un processus latent sous-jacent, https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/correlation-pluie05.gif tel que



A une date https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/correlation-pluie04.gif, on se demande ici s’il va pleuvoir dans la matinée, i.e. on veut

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/correlation-pluie07.gifqui peut encore s’écrire

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/correlation-pluie08.gifou, ce qui sera pleut être plus utile ensuite

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/correlation-pluie09.gifAvec notre modèle latent, cette probabilité peut s’écrire

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/correlation-pluie10.gifc’est à dire en notant https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/correlation-pluie13.gif la fonction de répartition (jointe) du vecteur aléatoire https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/correlation-pluie11.gif cette probabilité s’écrit
https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/correlation-pluie12.gifOr Wassily Hoeffding a montré qu’il était possible d’encadrer une fonction de répartition, avec les bornes suivantes: pour tout https://latex.codecogs.com/gif.latex?\boldsymbol{x}\in\mathbb{R}^d,

(en notant que ces bornes sont atteignables). Autrement dit, ici, la probabilité qu’il fasse beau 5 heures de suite est encadrée par la probabilité qu’il pleuve sur une heure… et 1.
Bref, on n’a pas grand chose à dire*.

Une possibilité est de supposer un modèle Gaussien pourhttps://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/correlation-pluie05.gif. Supposer les lois marginales n’est en rien contraignant. En revanche, supposer que la dynamique sous-jacente est Gaussienne est une hypothèse très forte. Et pour faire simple, supposons même que la dynamique soit auto-régressive (à l’ordre 1), i.e.
https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/correlation-pluie18.gifoù https://latex.codecogs.com/gif.latex?(\eta_t) est un bruit blanc, indépendant du processus latent. L’intérêt est que l’on a spécifié ainsi la structure de la matrice de variance covariance. En particulier, la matrice de corrélation du vecteur https://latex.codecogs.com/gif.latex?\boldsymbol{X}^\star est alors
On peut d’ailleurs visualiser facilement la forme de cette matrice en fonction de https://latex.codecogs.com/gif.latex?\rho,
Avec ce formalisme, on peut alors en déduire https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/correlation-pluie12.gif car on connaît très bien la fonction de répartition d’un vecteur gaussien multivarié via

On peut alors étudier l’impact de http://freakonometrics.blog.f ree.fr/public/perso3/slatex.gif (qui est associée à la probabilité d’avoir de la pluie pendant une heure en moyenne) et de la corrélation https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/correlation-pluie21.gif (qui est associée à la dynamique du processus). Pour simplifier, si on prend https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/slatex.gif tel que la probabilité de pluie soit de l’ordre de 25%, on a la probabilité d’avoir de la pluie suivante,

> s = .25
> r = .85
> i = 1:n
> I = matrix(rep(i,each=n),n,n)
> J = matrix(rep(i,n),n,n) 
> C = r^(abs(I-J))
> 1-pmnorm(rep(qnorm(1-s),n),mean=rep(0,n),varcov=C)[1]
[1] 0.458074

i.e. en faisant une boucle sur la valeur de la corrélation

Et plus généralement si on croise le seuil (i.e. la probabilité horaire de pleuvoir) et la probabilité, on a l’abaque suivante

On pourrait toutefois objecter que l’important n’est pas forcément qu’il ne pleuve pas de la matinée, mais qu’il ne pleuve pas quand ils sont dehors, i.e. (pour simplifier), à 8 heures, midi et 16 heures (c’est à dire toutes les 4 heures). On s’intéresse alors à
qui est toujours un vecteur Gaussien, mais avec cette fois une structure de corrélation un peu différente,
Dans ce cas, la probabilité de ne pas avoir de plus à ces trois moments de la journée est alors de la forme suivante

La corrélation joue alors un rôle relativement faible car il intervient à la puissance 4 dans la matrice d’auto-corrélation.
Dans ce cas, on notera que l’hypothèse d’indépendance n’est peut être pas absurde. Et avec une P.D.P. de 30%, sur toute la journée, on peut probablement se dire que la probabilité d’avoir de la pluie soit le matin, soit le midi, soit le soir sera de l’ordre de

> 1-(1-.3)^3
[1] 0.657

A tester !
* En fait si, on pourrait dire beaucoup de choses (et écrire une dizaine de billets) en particulier sur la borne supérieure, ici 1, dès lors que la probabilité horaire d’avoir de la pluie dépasse 20% – car la borne supérieure est


En fait, la borne inférieur n’est pas une fonction de répartition de manière générale. Toutefois, la borne est atteignable pour une certaine structure de dépendance. Ce qui explique que par la suite, avec un modèle Gaussien, on sera relativement loin de cette borne supérieure. La borne inférieure (i.e. la borne supérieure sur les fonction de répartition) est elle atteinte dès lors que les risques sont variables sont comonotones, c’est à dire dans le cas Gaussien si la corrélation entre deux instants (et entre tous les instants) vaut 1. La petite animation sur la matrice de corrélation permet de mieux comprendre la difficulté de comprendre le cas où la corrélation rho est négative: dans le cas où https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/rhofleche.gif1, on a une corrélation proche de 1 partout dans la matrice. En revanche, si https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/rhofleche.gif -1, on a des corrélations proches de -1 dans la moitié (environ) des cellules de la matrice, et +1 dans l’autre moitié. La borne supérieure sur la probabilité savoir de la pluie dans matinée correspond donc à une structure plus complexe que la borne inférieure, car on ne peut pas avoir une corrélation négative forte avec tous les lags.