Tag Archives: AR(1)

Forecasts with ARIMA Models

16/03/2016 Arthur Charpentier Leave a comment

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….

2014-E, Time Series

Variance of the Average of a Sequence

07/07/2014 Arthur Charpentier Leave a comment

In the case where $\{Y_1,\cdots,Y_n\}$ are i.i.d. random variables, then

$\text{Var}\left(\frac{1}{n}\sum_{t=1}^n Y_t\right)=\frac{\text{Var}(Y_t)}{n}$

Now, what if $\{Y_1,\cdots,Y_n\}$ are identically distributed, but no longer independent. What if we have an autoregressive process? Assume that

$Y_t=\phi Y_{t-1}+\varepsilon_t$

Then

$\text{Var}\left(\frac{1}{n}\sum_{t=1}^n Y_t\right) = \frac{1}{n^2}\left(\sum_{t=1}^n \text{Var}(Y_t)+\sum_{s\neq t}\text{Cov}(Y_s,Y_t) \right)$

can be written

$\text{Var}\left(\frac{1}{n}\sum_{t=1}^n Y_t\right) = \frac{1}{n^2}\left(n\gamma(0) + \sum_{h=1}^{n-1} 2 (n-h) \gamma(h)\right)$

Here, we will express the variance as a function of $\gamma(0)$ and $\phi$ , but it is possible to use also $\sigma^2$ , since, in the context of an $AR(1)$ ,

$\gamma(0)=\frac{\sigma^2}{1-\phi^2}$

Now, since $\gamma(h)=\phi^h \gamma(0)$ we get

$\text{Var}\left(\frac{1}{n}\sum_{t=1}^n Y_t\right) = \frac{\gamma(0)}{n^2}\left(n\gamma(0) + \sum_{h=1}^{n-1} 2 (n-h) \phi^h\right)$

which can be simplified, since

$\sum_{h=1}^{n-1} 2(n-h) \phi^h=2\phi^{n-1} \sum_{h=1}^{n-1} (n-h) \left(\frac 1\phi\right)^{n-h-1}=2\phi^{n-1} \frac {\partial}{\partial x}\left. \sum_{i=1}^{n} x^{n-i} \right\vert_{x=1/\phi}$

i.e.

$\sum_{h=1}^{n-1} 2(n-h) \phi^h=2\phi^{n-1} \frac{(n-1)\phi^{-(n+1)} - n\phi^{-n} + \phi^{-1}}{\phi^{-1}(\phi^{-1}-1)^2}=2 \frac{(n-1)\phi^{-1} - n + \phi^{n-1}}{(\phi^{-1}-1)^2}$

So, the variance of the mean can be writen as

$V=\mathrm{Var}\left(\frac{1}{n}\sum_{t=1}^n Y_t\right)=\frac{\gamma(0)}{n^2}\left[n + 2 \frac{(n-1)\phi^{-1} - n + \phi^{n-1}}{(\phi^{-1}-1)^2}\right]$

Observe that if $n$ is large enough,

$V=\frac{\gamma(0)}{n^2}\left[n + 2 \frac{(n-1)\phi^{-1} - n + \phi^{n-1}}{(\phi^{-1}-1)^2}\right]\sim \frac{\gamma(0)}{n}\frac{1+\phi}{1-\phi}$

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

$V=\mathrm{Var}\left(\frac{1}{n}\sum_{t=1}^n Y_t\right)=\frac{\gamma(0)}{n}\left[ \sum_{h=-n+1}^{n-1}\left(1-\frac{\vert h\vert}{n}\right)\rho(h) \right]$

or equivalently

$V=\mathrm{Var}\left(\frac{1}{n}\sum_{t=1}^n Y_t\right)=\frac{\gamma(0)}{n}\left[ 1+2\sum_{h=1}^{n-1}\left(1-\frac{h}{n}\right)\rho(h) \right]$

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

$V=\mathrm{Var}\left(\frac{1}{n}\sum_{t=1}^n Y_t\right)\sim\frac{\gamma(0)}{n}\left[ \sum_{h=-\infty}^{\infty}\rho(h) \right]$

In the context of an $AR(1)$ process, this can be writen

$\sum_{h=-\infty}^{\infty}\rho(h) = \frac{1+\phi}{1-\phi}$

Thus, we get the following well-known relationship

$V=\frac{\gamma(0)}{n}\cdot \frac{1+\phi}{1-\phi}$

In the case where $\{Y_1,\cdots,Y_n\}$ is an i.i.d. sequence, i.e. $\phi=0$ , then we get the relationship mentioned initially. And in the case of a random walk… unfortunately, we cannot use that relationship. But observe that

$V=\frac{1}{n}\text{Var}\left(\sum_{t=1}^n\sum_{h=1}^t\varepsilon_h\right)$

i.e.

$V=\frac{1}{n}\text{Var}\left(n\varepsilon_1+(n-1)\varepsilon_2+\cdots+\varepsilon_n\right)$ which can be written

$V=\frac{\text{Var}(\varepsilon)}{n}\sum_{h=1}^n h^2=\frac{(2n+1)(n+1)}{6n}\text{Var}(\varepsilon)$

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

MAT8181

Linear ‘Prediction’ for AR Time Series

25/03/2014 Arthur Charpentier Leave a comment

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,

$EL(Y\vert \boldsymbol{X})=\boldsymbol{\alpha}^{\star \text{\sffamily T}} \boldsymbol{X}$

where

$\boldsymbol{\alpha}^{\star}=\text{argmin}\left\{\mathbb{E}([Y-\boldsymbol{\alpha}^{\text{\sffamily T}} \boldsymbol{X}]^2)\right\}$

As discussed previously on this blog, we consider here on projection not on $\sigma\{\boldsymbol{X}\}$ (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 $\boldsymbol{X}$ is a (multivariate) Gaussian vector, $\mathcal{N}(\boldsymbol{0},\boldsymbol{\Sigma})$ , then

$\boldsymbol{\alpha}^{\star}=\boldsymbol{\Sigma}^{-1}\boldsymbol{\gamma}$

where $\boldsymbol{\gamma}$ is the vector $\text{cov}(Y,\boldsymbol{X})$ .

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

$X_{t}=\phi X_{t-1}+\varepsilon_t$

One observation (say $X_\tau$ ) is missing. We have here 3 questions:

what is the best linear prediction of $X_\tau$ given $X_{\tau-2}$ and $X_{\tau-1}$
what is the best linear prediction of $X_\tau$ given $X_{\tau+2}$ and $X_{\tau+1}$
what is the best linear prediction of $X_\tau$ given $X_{\tau-1}$ and $X_{\tau+1}$

Case 1. Here, we have to compute

$EL(X_{\tau}\vert X_{\tau-1}, X_{\tau-2})=\alpha_{-1}^\star X_{\tau-1}+\alpha_{-2}^\star X_{\tau-2}$

Since we have an AR(1) process, $\text{cov}(X_\tau,X_{\tau-1})=\phi$ , $\text{cov}(X_\tau,X_{\tau-2})=\phi^2$ and $\text{cov}(X_{\tau-1},X_{\tau-2})=\phi$ . Thus, from the relationship above

$\left(\begin{array}{ll}1 & \phi\\\phi & 1\end{array}\right)\left(\begin{array}{l}\alpha_{-2}^\star\\ \alpha_{-1}^\star\end{array}\right)=\left(\begin{array}{l}\phi^2\\\phi\end{array}\right)$

which can be written

$\left(\begin{array}{l}\alpha_{-2}^\star\\ \alpha_{-1}^\star\end{array}\right)=\left(\begin{array}{ll}1 & \phi\\\phi & 1\end{array}\right)^{-1}\left(\begin{array}{l}\phi^2\\\phi\end{array}\right)=\frac{1}{\phi^2-1}\left(\begin{array}{ll}-1 & \phi\\\phi & -1\end{array}\right)\left(\begin{array}{l}\phi^2\\\phi\end{array}\right)$

i.e. $(\alpha_{-2}^\star,\alpha_{-1}^\star)=(0,\phi)$ . Which makes sense actually: the AR(1) process is Markovian of order one, so

$EL(X_\tau\vert \underline{\boldsymbol{X}}_{\tau-1})=EL(X_\tau\vert X_{\tau-1})$

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

${}_{\tau}\widehat{X}_{\tau+1}=EL(X_\tau\vert \underline{\boldsymbol{X}}_{\tau-1})=\phi X_{\tau-1}$

So far, so good.

Case 2. Here, we have to compute

$EL(X_{\tau}\vert X_{\tau+1}, X_{\tau+2})=\alpha_{+1}^\star X_{\tau+1}+\alpha_{+2}^\star X_{\tau+2}$

Since we have an AR(1) process, $\text{cov}(X_\tau,X_{\tau+1})=\phi$ , $\text{cov}(X_\tau,X_{\tau+2})=\phi^2$ and $\text{cov}(X_{\tau+1},X_{\tau+2})=\phi$ . Thus, from the relationship above

$\left(\begin{array}{ll}1 & \phi\\\phi & 1\end{array}\right)\left(\begin{array}{l}\alpha_{+1}^\star\\ \alpha_{+2}^\star\end{array}\right)=\left(\begin{array}{l}\phi\\\phi^2\end{array}\right)$

i.e. $(\alpha_{+1}^\star,\alpha_{+2}^\star)=(\phi,0)$ .

Case 3. Finally, we have to compute

$EL(X_{\tau}\vert X_{\tau-1}, X_{\tau+1})=\alpha_{-1}^\star X_{\tau-1}+\alpha_{+1}^\star X_{\tau+1}$

One more time, since we have an AR(1) process, $\text{cov}(X_\tau,X_{\tau-1})=\phi$ , $\text{cov}(X_{\tau-1},X_{\tau+1})=\phi^2$ and $\text{cov}(X_{\tau},X_{\tau+1})=\phi$ . So here, the relationship above becomes

$\left(\begin{array}{ll}1 & \phi^2\\\phi^2 & 1\end{array}\right)\left(\begin{array}{l}\alpha_{-1}^\star\\ \alpha_{+1}^\star\end{array}\right)=\left(\begin{array}{l}\phi\\\phi\end{array}\right)$ Here, we can write

$\left(\begin{array}{l}\alpha_{-1}^\star\\ \alpha_{+1}^\star\end{array}\right)=\left(\begin{array}{ll}1 & \phi^2\\\phi^2 & 1\end{array}\right)^{-1}\left(\begin{array}{l}\phi\\\phi\end{array}\right)=\frac{1}{\phi^4-1}\left(\begin{array}{ll}-1 & \phi^2\\\phi^2 & -1\end{array}\right)\left(\begin{array}{l}\phi\\\phi\end{array}\right)$

i.e.

$\alpha_{-1}^\star=\frac{\phi}{\phi^2+1}=\alpha_{+1}^\star$

So finally, what we got here is

$EL(X_{\tau}\vert X_{\tau-2}, X_{\tau-1})=\phi X_{\tau-1}$

$EL(X_{\tau}\vert X_{\tau+1}, X_{\tau+2})=\phi X_{\tau+1}$

and

$EL(X_{\tau}\vert X_{\tau-1}, X_{\tau+1})=\frac{\phi}{\phi^2+1}\left(X_{\tau-1}+X_{\tau+1}\right)$

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

$\mathbb{E}\left([X_\tau-EL(X_{\tau}\vert\boldsymbol{X})]^2\right)$

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

Course, MAT8181, Time Series

Causal Autoregressive Time Series

21/01/2014 Arthur Charpentier 2 Comments

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)$ ,

$https://latex.codecogs.com/gif.latex?X_t=\phi%20X_{t-1}+\varepsilon_t$

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

$https://latex.codecogs.com/gif.latex?X_t%20=%20\sum_{h=0}^{+\infty}%20\phi^h%20\varepsilon_{t-h}$

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

$https://latex.codecogs.com/gif.latex?X_t%20=%20\sum_{h=0}^{+\infty}%20\phi^h%20\varepsilon_{t-h}$

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

$https://latex.codecogs.com/gif.latex?X_t%20=%20\sum_{h={\color{Red}%201}}^{+\infty}%20\frac{-1}{\phi^h}%20\varepsilon_{t{\color{Red}%20+h}}$

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 
0.908 

> 1/phi
[1] 0.9090909

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

$https://latex.codecogs.com/gif.latex?X_t=\frac{1}{\phi}X_{t-1}+\eta_t$

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

Course, MAT8181, Time Series

Visualizing Autoregressive Time Series

21/01/2014 Arthur Charpentier Leave a comment

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

e.g.

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

2012-H, ACT6420, Time Series

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

10/09/2012 Arthur Charpentier 1 Comment

Consider a time series, generated using

set.seed(1)
E=rnorm(240)
X=rep(NA,240)
rho=0.8
X[1]=0
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
u=241:(240+H)
P=X[240]*rhoest^(1:H)
s=sqrt(1/(sum((rhoest^(2*(1:300))))))*sd(X)

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

 
plot(1:240,X,xlab="",xlim=c(0,240+H),
ylim=c(-9.25,9),ylab="",type="l")
V=cumsum(rhoest^(2*(1:H)))*s
polygon(c(u,rev(u)),c(P+1.96*sqrt(V),
rev(P-1.96*sqrt(V))),col="yellow",border=NA)
polygon(c(u,rev(u)),c(P+1.64*sqrt(V),
rev(P-1.64*sqrt(V))),col="orange",border=NA)
lines(u,P,col="red")

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.

2012-H, Time Series

Simulation de séries temporelles

06/09/2012 Arthur Charpentier 1 Comment

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[1]=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[1]=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]}

Climate, Statistics

Trop de probabilités tue les probabilités

04/06/2011 Arthur Charpentier Leave a comment

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 la variable indicatrice indiquant s’il pleut, ou pas (1 s’il pleut, et 0 sinon) à l’instant . Supposons qu’il existe un processus latent sous-jacent, tel que

A une date , on se demande ici s’il va pleuvoir dans la matinée, i.e. on veut

qui peut encore s’écrire

ou, ce qui sera pleut être plus utile ensuite

Avec notre modèle latent, cette probabilité peut s’écrire

c’est à dire en notant la fonction de répartition (jointe) du vecteur aléatoire cette probabilité s’écrit
Or 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$ ,

où

et

(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 pour. 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.
où $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 car on connaît très bien la fonction de répartition d’un vecteur gaussien multivarié via

> library(mnormt)

On peut alors étudier l’impact de (qui est associée à la probabilité d’avoir de la pluie pendant une heure en moyenne) et de la corrélation (qui est associée à la dynamique du processus). Pour simplifier, si on prend 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

http://freakonometrics.hypotheses.org/files/2015/12/prooooobacommentaire.gif

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ù 1, on a une corrélation proche de 1 partout dans la matrice. En revanche, si -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.

Freakonometrics

Tag Archives: AR(1)

Forecasts with ARIMA Models

Variance of the Average of a Sequence

Linear ‘Prediction’ for AR Time Series

Causal Autoregressive Time Series

Visualizing Autoregressive Time Series

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

Simulation de séries temporelles

Trop de probabilités tue les probabilités

"sendo l'intento mio scrivere cosa utile a chi la intende…"