# Tents, Tweets, and Events: Ongoing Protests and Social Media

Our paper, entitled Tents, Tweets, and Events: The Interplay Between Ongoing Protests and Social Media, written with Marco Toledo Bastos (aka ) and Dan Mercea (aka) just appeared in the Journal of Communication

Recent protest movements have fuelled deliberations about the extent to which social media ignite protests. In this paper we compare time-series data of Twitter, Facebook, and onsite protest activity to test the hypothesis of Granger-causality between social media streams and protestors attending demonstrations during the Indignados in Spain, the Occupy movement in the U.S., and the Vinegar protests in Brazil. After applying a Gaussianization procedure to the time series, we confirmed the hypothesis that contentious communication on Twitter and Facebook was Granger-causal of onsite protest activity during the Indignados and the Occupy protests, with bidirectional causality between online and onsite protest activity in the Occupy series. The Vinegar protests in Brazil presented Granger-causality only between Facebook and Twitter and between protestors and injured or arrested protestors. The results indicate that the causal relationship between online and onsite political varies considerably across different socioeconomic contexts with different levels of Internet penetration.

# Variance of the Average of a Sequence

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

Pour le second devoir du cours ACT6420 (modèles de prévisions), le but est de prévoir des recherches sur Google, via https://www.google.com/trends/. Soit vous avez un mot clé qui vous intéresse, soit vous allez chercher une série extraite dans /ACT6420-TS/. Dans les bases qui ont été mises en ligne, il y a le mot clé gym, par exemple.

> report=read.table(
+ "http://freakonometrics.free.fr/ACT6420-TS/report-GYM.csv",
+ skip=4,header=TRUE,sep=",",nrows=548)

On va nettoyer un peu la base (en particulier les valeurs manquantes de la fin)

> tail(report)
Semaine gym
543 2014-05-25 - 2014-05-31  80
544 2014-06-01 - 2014-06-07  80
545 2014-06-08 - 2014-06-14  78
546 2014-06-15 - 2014-06-21  NA
547 2014-06-22 - 2014-06-28  NA
548 2014-06-29 - 2014-07-05  NA
> report=report[!is.na(report[,2]),]
> tail(report)
Semaine gym
540 2014-05-04 - 2014-05-10  79
541 2014-05-11 - 2014-05-17  80
542 2014-05-18 - 2014-05-24  79
543 2014-05-25 - 2014-05-31  80
544 2014-06-01 - 2014-06-07  80
545 2014-06-08 - 2014-06-14  78

Les données sont ici hebdomadaires, comme on peut le voir graphiquement

> hebdo=ts(report$gym,start=2004,frequency=52) > hebdo Time Series: Start = c(2004, 1) End = c(2014, 25) Frequency = 52 [1] 68 60 60 53 55 50 49 49 51 48 48 45 45 [14] 47 42 48 46 47 46 47 47 46 47 46 45 46 [27] 46 50 48 48 52 51 57 55 53 56 55 50 48 [40] 50 46 49 46 48 49 48 46 50 47 46 43 54 [53] 69 64 63 60 57 57 53 54 55 54 50 53 54 [66] 46 50 50 49 49 49 47 49 48 49 50 49 49 [79] 47 47 50 49 52 51 55 55 55 54 54 52 53 [92] 54 53 52 51 51 50 51 48 52 50 47 45 56 [105] 76 72 66 64 63 59 53 56 57 58 54 55 54 [118] 53 52 52 50 53 50 51 49 51 51 50 50 48 [131] 48 53 52 56 58 57 60 62 62 63 59 58 58 [144] 56 54 54 53 53 54 53 50 55 54 53 51 56 [157] 77 73 68 68 67 66 64 67 64 63 63 63 62 [170] 62 61 62 61 63 62 63 63 63 63 59 58 59 [183] 61 60 60 58 61 61 62 62 64 68 66 61 58 [196] 58 55 54 51 55 54 55 53 55 53 52 50 55 [209] 76 77 68 67 64 64 58 61 59 59 57 55 57 [222] 59 57 54 56 55 54 52 52 53 54 53 55 55 [235] 55 57 54 56 55 58 65 63 64 67 66 63 62 [248] 60 60 57 55 56 56 58 58 53 56 55 54 52 [261] 69 77 71 68 66 63 60 60 62 59 59 57 57 [274] 60 58 58 56 54 58 57 56 57 57 57 57 54 [287] 54 55 57 57 56 64 60 59 62 62 64 59 58 [300] 57 57 54 53 52 53 53 55 53 55 53 50 49 [313] 63 76 73 68 65 66 60 61 61 60 58 58 61 [326] 61 62 57 57 58 55 58 57 58 59 57 55 55 [339] 57 57 58 59 59 60 60 63 63 63 66 65 62 [352] 60 59 58 57 56 58 59 56 54 57 54 54 53 [365] 66 87 77 74 72 69 68 64 65 65 68 63 65 [378] 65 65 62 61 62 62 63 63 61 65 63 64 63 [391] 61 62 62 61 62 65 63 67 67 71 74 71 70 [404] 68 68 65 66 64 65 68 64 64 65 62 61 60 [417] 69 91 88 83 81 78 75 71 73 74 73 70 69 [430] 66 68 69 66 68 68 65 69 66 69 70 69 70 [443] 69 72 72 71 69 76 74 72 77 77 82 78 72 [456] 72 69 68 67 67 64 66 66 63 65 64 62 61 [469] 67 88 90 83 81 83 77 76 76 74 75 74 74 [482] 77 77 77 73 72 76 72 71 72 74 72 74 73 [495] 73 76 73 73 71 76 76 79 79 83 83 81 78 [508] 78 76 78 80 74 73 75 74 75 72 71 70 69 [521] 73 92 100 94 93 91 86 84 84 85 85 83 82 [534] 83 83 78 79 80 80 79 80 79 80 80 78 > plot(hebdo) Pour avoir des modèles plus simples, on peut agréger les données, pour les rendre mensuelles (par interpolation linéaire). La fonction à utiliser est ici  H2M=function(BASE){ X=BASE[,2] Y=BASE[,1] date1=substr(as.character(Y),1,10) date2=substr(as.character(Y),14,23) D1=as.Date(date1,"%Y-%m-%d") D2=as.Date(date2,"%Y-%m-%d") vm=vy=N=NA for(t in 1:length(D1)){ mois=seq(D1[t],D2[t],length=7) vm=c(vm,as.POSIXlt(mois)$mon+1)
vy=c(vy,as.POSIXlt(mois)$year+1900) N=c(N,rep(X[t],7))} N=N[-1]; vm=vm[-1]; vy=vy[-1] YM=vy*100+vm Z=tapply(N,as.factor(YM),mean) Zts=ts(as.numeric(Z),start=c(2004,1),frequency=12) return(Zts)} Si on utilise cette fonction sur nos données, on obtient > mensuel=H2M(report) > mensuel Jan Feb Mar Apr May Jun 2004 60.25000 50.75862 47.51613 45.66667 46.67742 46.00000 2005 63.22581 55.10714 53.03226 49.10000 48.45161 49.10000 2006 68.87097 57.10714 55.51613 51.86667 50.70968 49.93333 2007 70.74194 65.57143 62.87097 61.60000 62.77419 59.96667 2008 70.45161 60.79310 57.19355 56.13333 52.96774 54.30000 2009 70.35484 61.25000 58.19355 57.13333 56.80645 56.00000 2010 69.87097 61.78571 59.45161 58.76667 57.16129 56.40000 2011 76.58065 66.21429 65.22581 62.66667 62.51613 63.16667 2012 85.00000 73.82759 69.93548 67.76667 67.25806 69.46667 2013 84.93548 76.39286 75.00000 74.80000 72.67742 73.13333 2014 94.29032 84.96429 83.29032 79.80000 79.54839 79.00000 Jul Aug Sep Oct Nov Dec 2004 47.80645 53.67742 52.63333 47.77419 47.96667 47.61290 2005 48.41935 53.51613 53.43333 52.41935 50.20000 49.74194 2006 52.48387 59.77419 59.66667 54.12903 52.86667 54.35484 2007 59.87097 62.03226 63.10000 54.45161 54.30000 54.09677 2008 55.45161 62.16129 63.96667 57.41935 56.23333 56.09677 2009 55.96774 61.12903 60.16667 54.29032 53.60000 53.35484 2010 58.12903 61.64516 63.43333 57.67742 56.73333 56.48387 2011 61.80645 66.22581 70.86667 65.77419 65.23333 63.19355 2012 71.48387 75.06452 75.80000 67.22581 64.90000 65.12903 2013 73.51613 78.93548 79.73333 76.41935 73.93333 72.80645 2014 > ts.plot(mensuel) Cette base est maintenant utilisable pour le devoir. Le but est ici de faire de la prévision pour les 2 prochaines années, avec un intervalle de confiance. Mais j’en reparlerais par courriel. # Prévision de séries chronologiques Dans la seconde partie du cours de modèles de prévision, on quittera (un peu) les données individuelles pour parler de données chronologiques. Les slides de cette semaine (et probablement la semaine prochaine) sont en ligne. J’ai mis en ligne, en parallèle, quelques notes de cours sur les séries temporelles, qui pourront peut être servir de complément. Comme le disait Doug MartinTime series is the worst subject to teach. First, you have to teach the standard theory. Then, if you are beging honest, you have to tell the students ‘none of this stuff works, and this is what people really do'” (cité la semaine passée par David J. Thomson, au congrès de la SSC). On essayera de garder ça en mémoire tout au long du cours ! # Examen, Séries Chronologiques 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 ! # Stationarity of ARCH processes In the context of AR(1) processes, we spent some time to explain what happens when $\phi$ is close to 1. • if $\phi<1$ the process is stationary, • if $\phi=1$ the process is a random walk • if $\phi>1$ the process will explode Again, random walks are extremely interesting processes, with puzzling properties. For instance, $\text{Var}(X_{t+h}\vert X_t)\sim \sigma^2 h\rightarrow\infty$ as $h\rightarrow\infty$, and the process will cross the x-axis an infinite number of times… Recently, in the MAT8181 course, we studied carefully properties of the ARCH(1) process, especially when $\alpha\sim 1$. And again, what we get might be puzzling. Consider some ARCH(1) process $(\varepsilon_t)$, with a Gaussian noise, i.e. $\varepsilon_t=\sigma_t\cdot \eta_t$ where $\sigma^2_t=\omega+\alpha \varepsilon_{t-1}^2$ and $(\eta_t)$ is a sequence of i.i.d. $\mathcal{N}(0,1)$ variables. Here both $\omega$ and $\alpha$ have to be positive. Recall that $\mathbb{E}(\varepsilon_t)=0$ since $\mathbb{E}(\eta_t)=0$. Further $\text{var}(\varepsilon_t)=\mathbb{E}(\varepsilon_t^2)=\omega+\alpha \mathbb{E}(\varepsilon_{t-1}^2)$ since $\mathbb{E}(\eta_t^2)=1$, so the variance exists, and is constant only if $\alpha\in(0,1)$, and in that case $\sigma^2=\text{var}(\varepsilon_t)=\frac{\omega}{1-\alpha}\in[0,+\infty)$ Further, if $3\alpha^2<1$, then the fourth moment can be obtained, $\mathbb{E}(\varepsilon_t^4)=\frac{3\omega^2}{1-\alpha^2}\frac{1-\alpha^2}{1-3\alpha^2}$ since$\mathbb{E}(\eta_t^4)=3$. Now, if we get back on the property obtained while studying the variance, what does that mean if $\alpha=1$, or $\alpha>1$ ? If we look at simulations, we can generate an ARCH(1) process with $\alpha=2$ for instance. > n=600 > a=2 > w=0.2 > set.seed(1) > eta=rnorm(n) > epsilon=rnorm(n) > sigma2=rep(w,n) > for(t in 2:n){ + sigma2[t]=w+a*epsilon[t-1]^2 + epsilon[t]=eta[t]*sqrt(sigma2[t]) + } > plot(epsilon,type="l") In order to understand what’s going on, we should keep in mind that, what we good is that $\alpha$ has to lie in $(0,1)$ to be able to compute the second moment of $(\varepsilon_t)$. But it is possible to have a stationary process with infinite variance. And actually, this is what we have here. Write $\sigma^2_t=\omega+\alpha \varepsilon_{t-1}^2 = \omega+[\alpha \eta_{t-1}^2] \sigma_{t-1}^2$ and them, iterate $\sigma^2_t= \omega+[\alpha \eta_{t-1}^2] \left( \omega+[\alpha \eta_{t-2}^2] \sigma_{t-2}^2\right)$ and iterate again, and again, and again… $\sigma^2_t= \underbrace{\omega\left[ 1+\sum_{i=1}^h[\alpha \eta_{t-1}^2]\cdots [\alpha \eta_{t-i}^2] \right]}_{\Sigma_t(h)} +[\alpha \eta_{t-1}^2]\cdots[\alpha \eta_{t-h-1}^2]\sigma_{t-h-1}^2$ where $\Sigma_t(h)=\sum_{i=1}^h\underbrace{[\alpha \eta_{t-1}^2]\cdots [\alpha \eta_{t-i}^2]}_{u_i}$ Here, we have a sum of positive terms, and we can use the so-called Cauchy rule: define $\lambda=\text{limsup}\{ u_n^{1/n}\}$ then, if $\lambda<1$, the series $\sum u_n$ converges. Here, $u_n^{1/n}=\left[[\alpha \eta_{t-1}^2]\cdots [\alpha \eta_{t-n}^2]\right]^{1/n}$ which can also be written $u_n^{1/n}=\exp\left[\frac{1}{n}\sum_{i=1}^n\log[\alpha \eta_{t-i}^2]\right]$ and from the law of large numbers, since we have here a sum of i.i.d. terms, $u_n^{1/n}\rightarrow\exp\left[\mathbb{E}(\log[\alpha \eta^2])\right]$ So, if $\exp\left[\mathbb{E}(\log[\alpha \eta^2])\right]<1$, then $\Sigma_t(h)$ will have a limit when $h$ goes to infinity. The condition above can be written $\gamma=\mathbb{E}(\log[\alpha \eta^2]<0$ which is called Lyapunov coefficient. The equation $\exp\left[\mathbb{E}(\log[\alpha \eta^2])\right]=\alpha\exp\left[-\mathbb{E}(\log[\eta^2])\right]<1$ is a condition on $\alpha$. In the case where $\eta\sim\mathcal{N}(0,1)$, the numerical value of this upper bound is 3.56. > 1/exp(mean(log(rnorm(1e7)^2))) [1] 3.562517 In that case ($\gamma<0$), the variance may be infinite, but the series is stationary. On the other hand, if $\gamma>0$, then $\varepsilon_t^2$ will go to infinity almost surely, as $t$ goes to infinity. But in order to observe this difference, we need a lot of observations. For instance, with $\alpha=0.8$, and $\alpha=1.2$, we can easily see a difference. I do not say that it’s easy to see that the distribution above has an infinite variance, but still. Actually, if we consider Hill’s plot on the series above, on the tails of positive $\varepsilon_t$‘s > library(evir) > hill(epsilon) or on the tails of negative $\varepsilon_t$‘s > hill(-epsilon) we can see that the tail index is (strictly) smaller than 2 (meaning that the moment of order 2 does not exist). Why is it puzzling? Maybe because here, $(\varepsilon_t)$ is not weakly stationary (in the $L^2$ sense), but it is strongly stationary. Which is not the usual way weak and strong are related. This might be why we will not call this strong stationarity, but strict. # Inference for ARCH processes Consider some ARCH($p$) process, say ARCH($1$), $\varepsilon_t=\sigma_t\cdot \eta_t$ where $\sigma^2_t=\omega+\alpha \varepsilon_{t-1}^2$ with a Gaussian (strong) white noise $(\eta_t)$. > 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 $(\varepsilon_t)$ a ARCH($1$), then $(\varepsilon_t^2)$ is an AR($1$) process. So a first idea is to consider a regression, as we did for Gaussian AR($1$) > 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 $(\eta_t)$ is a Gaussian (strong) white noise, $\varepsilon_t\vert \underline{\boldsymbol{\epsilon}}_{t-1} \sim \mathcal{N}(0,\sigma^2_t)$ The likelihood is then $\mathcal{L}=\prod_{t=1}^T \frac{1}{\sqrt{2\pi \sigma^2_t}}\exp\left(-\frac{\varepsilon_t^2}{2\sigma_t}\right)$ and the log-likelihood is $\log \mathcal{L}=-\frac{1}{2} \sum_{t=1}^T \log (\sigma_t)-\frac{1}{2}\sum_{t=1}^T \frac{\varepsilon_t^2}{\sigma_t}$ And a natural idea is to define $(\widehat{\omega},\widehat{\alpha})\in\text{argmax}\{\log \mathcal{L}(\omega,\alpha)\}$ 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 $\alpha$ 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  $\alpha$ 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($2$) time series

$\varepsilon_t=\sigma_t\cdot \eta_t$

where now

$\sigma^2_t=\omega+\alpha_1 \varepsilon_{t-1}^2+\alpha_2 \varepsilon_{t-2}^2$

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

$\log \mathcal{L}=-\frac{1}{2} \sum_{t=1}^T \log (\sigma_t)-\frac{1}{2}\sum_{t=1}^T \frac{\varepsilon_t^2}{\sigma_t}$

and we can define

$(\widehat{\omega},\widehat{\alpha}_1,\widehat{\alpha}_2)\in\text{argmax}\{\log \mathcal{L}(\omega,\alpha_1,\alpha_2)\}$

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($1$,$1$), $\varepsilon_t=\sigma_t\cdot \eta_t$ where now $\sigma^2_t=\omega+\alpha \varepsilon_{t-1}^2+\beta \sigma_{t-1}$ Again, maximum likelihood techniques can be used. Actually, we can also code Fisher-Scoring algorithm, since (in a very general context) $\frac{\partial \log \mathcal{L}}{\partial \boldsymbol{\theta}}=-\frac{1}{2}\sum_{t=1}^T \frac{\partial \sigma_t^2}{\partial \boldsymbol{\theta}}\cdot \frac{1}{\sigma_t^2}\left(1-\frac{\varepsilon^2_t}{\sigma_t^2}\right)$ with here $\boldsymbol{\theta}=(\omega,\alpha,\beta)$. 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 # 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 $(1-L^{12})\Phi(L)[X_t-\mu]=\Theta(L)\varepsilon_t$ If we forget about the autoregressive and the moving average component, we can estimate $(1-{\color{red}\alpha}L^{12})[X_t-\mu]=\varepsilon_t$ If there is a seasonal unit root then ${\color{red}\alpha}$ 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 $VAR(1)$ model on the quarterly temperature? Something like $\left(\begin{matrix}X_t^{(Q1)}\\ X_t^{(Q2)\\ X_t^{(Q3)\\ X_t^{(Q4)\end{matrix}\right)= A \left(\begin{matrix}X_{t-1}^{(Q1)}\\ X_{t-1}^{(Q2)\\ X_{t-1}^{(Q3)\\ X_{t-1}^{(Q4)\end{matrix}\right)+\left(\begin{matrix}\varepsilon_{t-1}^{(Q1)}\\ \varepsilon_{t-1}^{(Q2)\\ \varepsilon_{t-1}^{(Q3)\\ \varepsilon_{t-1}^{(Q4)\end{matrix}\right)$ i.e. $\boldsymbol{Y}_t=A \boldsymbol{Y}_{t-1}+\boldsymbol{\varepsilon}_t$ where $A$ is some matrix $4\times 4$. 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 $A$ 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 $PAR(1)$ model, discussed by Paap and Franses $\boldsymbol{\Phi}_0 \boldsymbol{Y}_t=\boldsymbol{\Phi}_1 \boldsymbol{Y}_{t-1}+\boldsymbol{\varepsilon}_t$ where $\boldsymbol{\Phi}_0=\left(\begin{array}{cccc}1&0&0&0\\ -\phi_{1,2}&1&0&0\\ 0&-\phi_{1,3}&1&0\\ 0&0&-\phi_{1,4}&1\end{array}\right)$ and $\boldsymbol{\Phi}_1=\left(\begin{array}{cccc}0&0&0&\phi_{1,1}\\ 0&0&0&0\\ 0&0&0&0\\ 0&0&0&0\end{array}\right)$ Keep in mind that this is a $VAR(1)$ model, since $\boldsymbol{Y}_t=\underbrace{\boldsymbol{\Phi}_0^{-1} \boldsymbol{\Phi}_1}_A \boldsymbol{Y}_{t-1}+\boldsymbol{\varepsilon}_t$ 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 $|\boldsymbol{\Phi}_0-\boldsymbol{\Phi}_1 z |=(1-\phi_{1,1}\phi_{1,2}\phi_{1,3}\phi_{1,4}\cdot z)$ so there is a (seasonal) unit root if $\phi_{1,1}\phi_{1,2}\phi_{1,3}\phi_{1,4}=1$ 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 $(1-L^4)$ has four root, in $\mathbb{C}$, $\{+1,+i,-1,-i\}$ since $(1-L^4)=(1-L)(1+L)(1+L^2)=(1-L)(1+L)(1-iL)(1+iL)$ If we get back to monthly data, $(1-L^{12})$ has twelve roots, $\left\{\pm 1, \ \pm i, \ \pm\frac{1\pm \sqrt{3}i}{2},\pm\frac{\sqrt{3}\pm i}{2}\right\}$ 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 $\left\{\frac{\pi}{6},\frac{\pi}{3},\frac{\pi}{2},\frac{2\pi}{3},\frac{5\pi}{6},\pi\right\}$ 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). # 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 $(X_t)$, 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_t=X_t-[a+bt]$

> 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 $Z_t=(1-L^{12})Y_t$, 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

$(1-\phi L){Z_t}=(1-\phi L)\underbrace{(1-L^{12})Y_t}=\varepsilon_t$

for some white noise $(\varepsilon_t)$. This suggest the following SARIMA structure on $(Y_t)$,

> 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

$(1-\phi L){(1-{\color{red}\phi_{12}}L^{12})Y_t}=\varepsilon_t$

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 $p$, then we get a model with $p$ 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 $(X_t)_{t\in\mathbb{Z}}$ be a stationary time series, and $(a_t)_{t\in\mathbb{Z}}$ a sequence of real numbers such that

$\sum_{t\in\mathbb{Z}} \vert a_t\vert < \infty$

then the time series $(Y_t)_{t\in\mathbb{Z}}$ defined as

$Y_t=\sum_{k\in\mathbb{Z}} a_k X_{t-k}= \left( \sum_{k\in\mathbb{Z}} a_k L^k\right) X_t$

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

$\gamma_Y(h)=\text{cov}(Y_t,Y_{t-k})=\sum_{r,s\in\mathbb{Z}} a_r a_s \underbrace{\gamma_X(h+r-s)}_{\text{cov}(X_{t-r},X_{t-s+h})}$

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 $f_Y(\omega)=\vert A(e^{i\omega})\vert^2f_X(\omega)$ where $A(\cdot)$ 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 $(Y_t)$ is (weakly) stationnary then $Y_{t}=\sum _{{j=0}}^{\infty }\psi_{j}\varepsilon _{{t-j}}+\xi _{t}$ where $(\varepsilon _{{t}})$ is the innovation process, and where $(\xi _{{t}})$ is some deterministic series (just to get a result as general as possible). Observe that $\sum _{{j=0}}^{{\infty }}|\psi_{{j}}|^{2} < \infty$ 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 $\Delta Y_{t}=(1-L) Y_t=\sum _{{j=0}}^{\infty }\psi_{j}\varepsilon _{{t-j}}+\xi=\Psi(L)\varepsilon _{{t}}+\xi$can be represented as a linear trend $+$ a random walk $+$ a stationary remaining term i.e. $Y_{t}=\underbrace{Y_0 + \xi t }+\underbrace{\Psi(1)\sum_{j=1}^t\varepsilon _{{i}}}+\underbrace{\tilde\Psi(L)\varepilon_0-\tilde\Psi(L)\varepsilon_t}$ where $\tilde\Psi(\cdot)$ is the polynomial with terms $\tilde\psi_j$, where $\tilde\psi_j =\sum_{i=j+1}^\infty\psi_i$ 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 $Y_t=\alpha+\beta t+\varphi Y_{t-1}+\varepsilon_t$ and we would like to test if $\varphi=1$ (or not). We can write the previous representation as $\Delta Y_t=\alpha+\beta t+[\varphi-1] Y_{t-1}+\varepsilon_t$ 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 $\varphi-1=0$. 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)

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

$\Delta Y_t=\alpha+\beta t+[\varphi-1] Y_{t-1}+\psi \Delta Y_{t-1}+\varepsilon_t$

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 $ARMA(1,1)$ process,

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

where $(\varepsilon_t)$ is some white noise, and assume further that $\theta+\phi\neq0$.

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

$X_t=\phi X_{t-1}+u_t$

The estimator of $\phi$ – 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

$u_t=\varepsilon_t+\theta\varepsilon_{t-1}$

meaning that $(u_t)$ should be some $MA(1)$ process. And

$\rho(1)=\frac{\theta}{1+\theta^2}$

i.e. $\theta$ 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,

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

$\gamma(1)=\phi\gamma(0)+\theta\sigma^2$

and

$\gamma(2)=\phi\gamma(1)$

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 $MA(q)$ 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 $\phi\in(-1,1)$, > 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
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 $AR(p)$ 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 $MA(2)$ time series

$X_t=\varepsilon_t +\theta_1 \varepsilon_{t-1}+\theta_2\varepsilon_{t-2}$

for some white noise $(\varepsilon_t)$.

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

$\rho(1)=\frac{\theta_1+\theta_1\theta_2}{1+\theta_1^2+\theta_2^2}$

$\rho(2)=\frac{\theta_2}{1+\theta_1^2+\theta_2^2}$

with $\rho(h)=0$ when $h\geq0$. We also have the following relationship on the variance of the process

$\text{Var}(X_t)=\sigma^2\cdot(1+\theta_1^2+\theta_2^2)$

With those three equations, for three unknown parameters, $\theta_1$$\theta_2$ and $\sigma$, 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 $(\varepsilon_t)$, we only observe our series $(X_t)$. But we can try to rebuild that series (call it $(u_t)$ 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 $u_t =X_t-\left(\theta_1u_{t-1}+\theta_2un_{t-2}\right)$ and then, we can use least square techniques $\min_{(\theta_1,\theta_2)}\left\{\sum_{i=1}^tu_i^2\right\}$ 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
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 $\boldsymbol{Y}=(Y_1,\cdots,Y_t)$ 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 $AR(p)$ time series.