Category Archives: Time Series

Testing for a causal effect (with 2 time series)

19/02/2020 Arthur Charpentier 7 Comments

A few days ago, I came back on a sentence I found (in a French newspaper), where someone was claiming that

“… an old variable explains 85% of the change in a new variable. So we can talk about causality”

and I tried to explain that it was just stupid : if we consider the regression of the temperature on day $t+1$ against the number of cyclist on day $t$ , the $R^2$ exceeds 80%… but it is hard to claim that the number of cyclists on specific day will actually cause the temperature on the next day…

Nevertheless, that was frustrating, and I was wondering if there was a clever way to test for causality in that case. A popular one is Granger causality (I can mention a paper we published a few years ago where we use such a test, Tents, Tweets, and Events: The Interplay Between Ongoing Protests and Social Media). To explain that test, consider a bivariate time series (just like the one we have here), $\boldsymbol{z}_t=(x_t,y_t)$ , and consider some bivariate autoregressive model
${\begin{bmatrix}x_{t}\\y_{t}\end{bmatrix}}={\begin{bmatrix}c_{1}\\c_{2}\end{bmatrix}}+{\begin{bmatrix}a_{1,1}&\textcolor{red}{a_{1,2}}\\\textcolor{blue}{a_{2,1}}&a_{2,2}\end{bmatrix}}{\begin{bmatrix}x_{t-1}\\y_{t-1}\end{bmatrix}}+{\begin{bmatrix}u_{t}\\v_{t}\end{bmatrix}}$ where $\boldsymbol{\varepsilon}_t=(u_t,v_t)$ is some bivariate white noise, in the sense that (i) $\mathbb{E} (\boldsymbol{\varepsilon}_{t})=\boldsymbol{0}$ (the noise is centered) (ii) $\mathbb{E} (\boldsymbol{\varepsilon}_{t}\boldsymbol{\varepsilon}_{t}^\top)=\Omega$ , so the variance matrix is constant, but possibly non-diagonal (iii) $\mathbb{E} (\boldsymbol{\varepsilon}_{t}\boldsymbol{\varepsilon}_{t-h}^\top)=\boldsymbol{0}$ for all $h\neq 0$ . Note that we can use the simplified expression ${\boldsymbol{z}_t=\boldsymbol{c}+\boldsymbol{A}\boldsymbol{z}_{t-1}+\boldsymbol{\varepsilon}_t}$ Now, Granger test is based on several quantities. With off-diagonal terms of matrix $\Omega$ , we have a so-called instantaneous causality, and since $\Omega$ is symmetry, we will write $x\leftrightarrow y$ . With off-diagonal terms of matrix $\boldsymbol{A}$ , we have a so-called lagged causality, with either $\textcolor{blue}{x\rightarrow y}$ or $\textcolor{red}{x\leftarrow y}$ (and possibly both, if both terms are significant).

So I wanted to try on my two-variable problem.

df = read.csv("cyclistsTempHKI.csv")
dfts = cbind(C=ts(df$cyclists,start = c(2014, 1,2), frequency = 365),
             T=ts(df$meanTemp,start = c(2014, 1,2), frequency = 365))
library(vars)

I now have “time series” objects, and we can fit a VAR model,

var2 = VAR(dfts, p = 1, type = "const")
coefficients(var2)
$C
         Estimate   Std. Error   t value      Pr(&gt;|t|)
C.l1    0.8684009   0.02889424 30.054460 8.080226e-107
T.l1   70.3042012  20.07247411  3.502518  5.102094e-04
const 807.6394001 187.75472482  4.301566  2.110412e-05
 
$T
           Estimate   Std. Error   t value     Pr(&gt;|t|)
C.l1   0.0003865391 6.257596e-05  6.177118 1.540467e-09
T.l1   0.6611135594 4.347074e-02 15.208241 6.086394e-42
const -1.6413074565 4.066184e-01 -4.036481 6.446018e-05

For instant, we can run a causality, to test if the number of cyclists can cause the temperature (on the next day)

causality(var2, cause = "C")
$Granger
 
	Granger causality H0: C do not Granger-cause T
 
data:  VAR object var2
F-Test = 38.157, df1 = 1, df2 = 842, p-value = 1.015e-09

Here, we should clearly reject $H_0$ , which is that there is no causal effect. Which is the way statisticians say that there should be some causal effect between the number of cyclist and the temperature…

So clearly, something is wrong here. Either it is some sort of superpower that cyclists are not aware of. Or this test that was used for forty years (Clive Granger even got a Nobel price for it) is not working. Or we missed something. Actually… I think we missed something here. Possibly because the series are not stationary. We can almost see it with

Phi = matrix(c(coefficients(var2)$C[1:2,1],coefficients(var2)$T[1:2,1]),2,2)
eigen(Phi)
eigen() decomposition
$values
[1] 0.9594810 0.5700335

where the highest eigenvalue is very close to one. But actually, we look here at the temperature…

plot(dfts)

so, at least, we should expect some seasonal unit root here. So let us use two techniques. The first one is a classical one-year difference, $\Delta_{365}\boldsymbol{z}_t=\boldsymbol{z}_t-\boldsymbol{z}_{t-365}$

var2 = VAR(diff(dfts,365), p = 1, type = "const")
coefficients(var2)
$C
          Estimate   Std. Error   t value     Pr(&gt;|t|)
C.l1     0.8376424   0.07259969 11.537823 1.993355e-16
T.l1    42.2638410  28.58783276  1.478386 1.449076e-01
const -507.5514795 219.40240747 -2.313336 2.440042e-02
 
$T
         Estimate   Std. Error   t value     Pr(&gt;|t|)
C.l1  0.000518209 0.0003277295 1.5812096 1.194623e-01
T.l1  0.598425288 0.1290511945 4.6371154 2.162476e-05
const 0.547828079 0.9904263469 0.5531235 5.823804e-01

The test on the fited VAR model yields

causality(var2, cause = "C") 
$Granger
 
	Granger causality H0: C do not Granger-cause T
 
data:  VAR object var2
F-Test = 2.5002, df1 = 1, df2 = 112, p-value = 0.1167

i.e., with a 11% $p$ -value, we should reject the assumption that the number of cyclists cause the temperature (on the next day), and actually, we should also reject the other way

causality(var2, cause = "T") 
$Granger
 
	Granger causality H0: T do not Granger-cause C
 
data:  VAR object var2
F-Test = 2.1856, df1 = 1, df2 = 112, p-value = 0.1421

Nevertheless, if we look at the instantaneous causality, this one makes more sense

$Instant
 
	H0: No instantaneous causality between: T and C
 
data:  VAR object var2
Chi-squared = 13.081, df = 1, p-value = 0.0002982

The second idea would be to use a one day difference, $\Delta_{1}\boldsymbol{z}_t=\boldsymbol{z}_t-\boldsymbol{z}_{t-1}$ and to fit a VAR model on that one

VARselect(diff(dfts,1), lag.max = 4, type="const")
$selection
AIC(n)  HQ(n)  SC(n) FPE(n) 
     3      3      2      3

but on that one, a VAR(1) model – with only one lag – might not be sufficient. It might be better to consider a VAR(3)

var2 = VAR(diff(dfts,1), p = 3, type = "const")

and on that one, one more time, we should reject the causal effect of the number of cyclists on the temperature (on the next day)

causality(var2, cause = "C")  
$Granger
 
	Granger causality H0: C do not Granger-cause T
 
data:  VAR object var2
F-Test = 0.67644, df1 = 3, df2 = 828, p-value = 0.5666

and this time, there could be a (lagged) causal effect of the temperature on the number of cyclists

causality(var2, cause = "T")  
$Granger
 
	Granger causality H0: T do not Granger-cause C
 
data:  VAR object var2
F-Test = 7.7981, df1 = 3, df2 = 828, p-value = 3.879e-05
 
$Instant
 
	H0: No instantaneous causality between: T and C
 
data:  VAR object var2
Chi-squared = 55.83, df = 1, p-value = 7.905e-14

but nothing instantaneously… So it looks like Granger causality performs well on that one !

Publications, Time Series

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

03/03/2015 Arthur Charpentier 1 Comment

Our paper, entitled Tents, Tweets, and Events: The Interplay Between Ongoing Protests and Social Media, written with Marco Toledo Bastos (aka @toledobastos) and Dan Mercea (aka @DanMercea) 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.

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

2014-E, ACT6420, Time Series

Google, et prévision

10/06/2014 Arthur Charpentier Leave a comment

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.

2014-E, ACT6420, Course, Time Series

Prévision de séries chronologiques

09/06/2014 Arthur Charpentier Leave a comment

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 Martin “Time 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 !

MAT8181, Time Series

Examen, Séries Chronologiques

22/04/2014 Arthur Charpentier Leave a comment

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 !

MAT8181, Time Series

Stationarity of ARCH processes

06/04/2014 Arthur Charpentier Leave a comment

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.

MAT8181, Time Series

Inference for ARCH processes

02/04/2014 Arthur Charpentier 3 Comments

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

MAT8181, Time Series

Seasonal Unit Roots

26/03/2014 Arthur Charpentier 2 Comments

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)

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

Seasonal, or periodic, time series

20/03/2014 Arthur Charpentier 2 Comments

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 $https://latex.codecogs.com/gif.latex?p$ , then we get a model with $https://latex.codecogs.com/gif.latex?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.

MAT8181, Time Series

Filtering a Stationary Time Series

13/03/2014 Arthur Charpentier 1 Comment

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…

MAT8181, Tests, Time Series

Identification of ARMA processes

19/02/2014 Arthur Charpentier 3 Comments

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

MAT8181, Time Series

Temperatures Series as Random Walks

12/02/2014 Arthur Charpentier 3 Comments

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

MAT8181, Time Series

Unit Root Tests

12/02/2014 Arthur Charpentier 5 Comments

This week, in the MAT8181 Time Series course, we’ve discussed unit root tests. According to Wold’s theorem, if $https://latex.codecogs.com/gif.latex?(Y_t)$ is (weakly) stationnary then

$https://latex.codecogs.com/gif.latex?%20%20%20%20%20Y_{t}=\sum%20_{{j=0}}^{\infty%20}\psi_{j}\varepsilon%20_{{t-j}}+\xi%20_{t}$

where $https://latex.codecogs.com/gif.latex?%20(\varepsilon%20_{{t}})$ is the innovation process, and where $https://latex.codecogs.com/gif.latex?%20(\xi%20_{{t}})$ is some deterministic series (just to get a result as general as possible). Observe that

$https://latex.codecogs.com/gif.latex?\sum%20_{{j=0}}^{{\infty%20}}|\psi_{{j}}|^{2}%20%3C%20\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

$https://latex.codecogs.com/gif.latex?%20%20%20%20%20\Delta%20Y_{t}=(1-L)%20Y_t=\sum%20_{{j=0}}^{\infty%20}\psi_{j}\varepsilon%20_{{t-j}}+\xi=\Psi(L)\varepsilon%20_{{t}}+\xi$ can be represented as

a linear trend $https://latex.codecogs.com/gif.latex?+$ a random walk $https://latex.codecogs.com/gif.latex?+$ a stationary remaining term

i.e.

$https://latex.codecogs.com/gif.latex?%20Y_{t}=\underbrace{Y_0%20+%20\xi%20t%20%C2%A0}+\underbrace{\Psi(1)\sum_{j=1}^t\varepsilon%20_{{i}}}+\underbrace{\tilde\Psi(L)\varepilon_0-\tilde\Psi(L)\varepsilon_t}$

where $https://latex.codecogs.com/gif.latex?%20\tilde\Psi(\cdot)$ is the polynomial with terms $https://latex.codecogs.com/gif.latex?%20\tilde\psi_j$ , where

$https://latex.codecogs.com/gif.latex?%20\tilde\psi_j%20=\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

$https://latex.codecogs.com/gif.latex?%20Y_t=\alpha+\beta%20t+\varphi%20Y_{t-1}+\varepsilon_t$

and we would like to test if $https://latex.codecogs.com/gif.latex?%20\varphi=1$ (or not). We can write the previous representation as

$https://latex.codecogs.com/gif.latex?%20\Delta%20Y_t=\alpha+\beta%20t+[\varphi-1]%20Y_{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 $https://latex.codecogs.com/gif.latex?%20\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)
> adf.test(X,k=0)

	Augmented Dickey-Fuller Test

data:  X
Dickey-Fuller = -2.0433, Lag order = 0, p-value = 0.5576
alternative hypothesis: stationary

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

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

(but I think it is a coincidence).

Augmented Dickey Fuller

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

$https://latex.codecogs.com/gif.latex?%20\Delta%20Y_t=\alpha+\beta%20t+[\varphi-1]%20Y_{t-1}+\psi%20\Delta%20Y_{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.

MAT8181, Time Series

Inference for ARMA(p,q) Time Series

30/01/2014 Arthur Charpentier 1 Comment

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

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

where $https://latex.codecogs.com/gif.latex?(\varepsilon_t)$ is some white noise, and assume further that $https://latex.codecogs.com/gif.latex?\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,

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

The estimator of $https://latex.codecogs.com/gif.latex?\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

$https://latex.codecogs.com/gif.latex?u_t=\varepsilon_t+\theta\varepsilon_{t-1}$

meaning that $https://latex.codecogs.com/gif.latex?(u_t)$ should be some $https://latex.codecogs.com/gif.latex?MA(1)$ process. And

$https://latex.codecogs.com/gif.latex?\rho(1)=\frac{\theta}{1+\theta^2}$

i.e. $https://latex.codecogs.com/gif.latex?\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,

$https://latex.codecogs.com/gif.latex?\gamma(0)=\frac{1+\theta^2+2\phi\theta}{1-\phi^2}\cdot\sigma^2$

$https://latex.codecogs.com/gif.latex?\gamma(1)=\phi\gamma(0)+\theta\sigma^2$

and

$https://latex.codecogs.com/gif.latex?\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 $https://latex.codecogs.com/gif.latex?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 $https://latex.codecogs.com/gif.latex?\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
function gradient 
     118       NA 

$convergence
[1] 0

$message
NULL

$outer.iterations
[1] 2

$barrier.value
[1] -1.536358e-05

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

Freakonometrics

Category Archives: Time Series

Testing for a causal effect (with 2 time series)

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

Variance of the Average of a Sequence

Google, et prévision

Prévision de séries chronologiques

Examen, Séries Chronologiques

Stationarity of ARCH processes

Inference for ARCH processes

Seasonal Unit Roots

Seasonal, or periodic, time series

Filtering a Stationary Time Series

Identification of ARMA processes

Temperatures Series as Random Walks

Unit Root Tests

Inference for ARMA(p,q) Time Series

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