Kendall’s function for copulas

As mentioned in the course on copulas, a nice tool to describe dependence it Kendall’s cumulative function. Given a random pair http://freakonometrics.hypotheses.org/files/2015/12/conc-19.gif with distribution  http://freakonometrics.hypotheses.org/files/2015/12/conc-17.gif, define random variable http://freakonometrics.hypotheses.org/files/2015/12/conc-30.gif. Then Kendall’s cumulative function is

http://freakonometrics.hypotheses.org/files/2015/12/kendall-01.gif

Genest and Rivest (1993) introduced that function to choose among Archimedean copulas (we’ll get back to this point below).

From a computational point of view, computing such a function can be done as follows,

  • for all http://freakonometrics.hypotheses.org/files/2015/12/kendall-02.gif, compute http://freakonometrics.hypotheses.org/files/2015/12/kendall-03.gif as the proportion of observation in the lower quadrant, with upper corner http://freakonometrics.hypotheses.org/files/2015/12/kendall-4.gif, i.e.

http://freakonometrics.hypotheses.org/files/2015/12/kendall-06.gif

  • then compute the cumulative distribution function of http://freakonometrics.hypotheses.org/files/2015/12/kendall-03.gif‘s.

To visualize the construction of that cumulative distribution function, consider the following animation

Thus, here the code to compute simply that cumulative distribution function is

n=nrow(X)
i=rep(1:n,each=n)
j=rep(1:n,n)
S=((X[i,1]>X[j,1])&(X[i,2]>X[j,2]))
Z=tapply(S,i,sum)/(n-1)

The graph can be obtain either using

plot(ecdf(Z))

or

plot(sort(Z),(1:n)/n,type="s",col="red")

The interesting point is that for an Archimedean copula with generator http://freakonometrics.hypotheses.org/files/2015/12/kendall-7.gif, then Kendall’s function is simply

http://freakonometrics.hypotheses.org/files/2015/12/kendall-8.gifIf we’re too lazy to do the maths, at least, it is possible to compute those functions numerically. For instance, for Clayton copula,

h=.001
phi=function(t){(t^(-alpha)-1)}
dphi=function(t){(phi(t+h)-phi(t-h))/2/h}
k=function(t){t-phi(t)/dphi(t)}
Kc=Vectorize(k)

Similarly, let us consider Gumbel copula,

phi=function(t){(-log(t))^(theta)}
dphi=function(t){(phi(t+h)-phi(t-h))/2/h}
k=function(t){t-phi(t)/dphi(t)}
Kg=Vectorize(k)

If we plot the empirical Kendall’s function (obtained from the sample), with different theoretical ones, derived from Clayton copulas (on the left, in blue) or Gumbel copula (on the right, in purple), we have the following,

http://freakonometrics.hypotheses.org/files/2015/12/kendall-function-anim.gif

Note that the different curves were obtained when Clayton copula has Kendall’s tau equal to 0, .1, .2, .3, …, .9, 1, and similarly for Gumbel copula (so that Figures can be compared). The following table gives a correspondence, from Kendall’s tau to the underlying parameter of a copula (for different families)

as well as Spearman’s rho,


To conclude, observe that there are two important particular cases that can be identified here: the case of perfect dependent, on the first diagonal when http://freakonometrics.hypotheses.org/files/2015/12/kennnn-04.gif, and the case of independence, the upper green curve, http://freakonometrics.hypotheses.org/files/2016/10/kennnnn-05.gif. It should also be mentioned that it is also common to plot not function http://freakonometrics.hypotheses.org/files/2015/12/kennnn-01.gif, but function http://freakonometrics.hypotheses.org/files/2015/12/kennnn-02.gif, defined as http://freakonometrics.hypotheses.org/files/2015/12/kennnn-03.gif,

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

Consider a time series, generated using

set.seed(1)
E=rnorm(240)
X=rep(NA,240)
rho=0.8
X[1]=0
for(t in 2:240){X[t]=rho*X[t-1]+E[t]}

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

rhoest=0.9; H=260
u=241:(240+H)
P=X[240]*rhoest^(1:H)
s=sqrt(1/(sum((rhoest^(2*(1:300))))))*sd(X)

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

 
plot(1:240,X,xlab="",xlim=c(0,240+H),
ylim=c(-9.25,9),ylab="",type="l")
V=cumsum(rhoest^(2*(1:H)))*s
polygon(c(u,rev(u)),c(P+1.96*sqrt(V),
rev(P-1.96*sqrt(V))),col="yellow",border=NA)
polygon(c(u,rev(u)),c(P+1.64*sqrt(V),
rev(P-1.64*sqrt(V))),col="orange",border=NA)
lines(u,P,col="red")
Here, forecasts can be derived, with any kind of possible autoregressive coefficient, from 0.7 to 1. I.e. we can chose to model the time series either with a stationary, or an integrated series,

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

Les tests de non-stationnarité (racine unité)

Pour commencer la modélisation d’une série temporelle, la première étape est de savoir si on peut la considérer comme stationnaire. L’alternative étant qu’elle soit intégrée. Le test le plus classique est le test de Dickey-Fuller.De manière générale, considérons un modèle autorégressive à l’ordre 1,

http://freakonometrics.blog.free.fr/public/perso6/adf-1.gif

On peut éventuellement réécrire ce processus sous la forme

http://freakonometrics.blog.free.fr/public/perso6/acf-18.gif

http://freakonometrics.blog.free.fr/public/perso6/adf-21.gif. On aura un processus non-stationnaire dès lors que http://freakonometrics.blog.free.fr/public/perso6/adf-20.gif ou encore http://freakonometrics.blog.free.fr/public/perso6/adf-16.gif. Si on n’a pas de tendance linéaire, tester http://freakonometrics.blog.free.fr/public/perso6/acf-22.gif ou http://freakonometrics.blog.free.fr/public/perso6/adf-26.gif peut se faire avec un test de Student. L’idée de Dickey-Fuller est de généraliser le test de Student, en rajoutant cette tendance linéaire (pour commencer).

Afin d’illustrer l’utilisation de ce test, commençons par un cas “simple”: on va intégrer un bruit blanc simulé.

>  set.seed(1)
>  E=rnorm(240)
>  X=cumsum(E)
>  plot(X,type="l")

On peut regarder s’il y a une constante non nulle (on parlera de modèle avec drift) voire une tendance linéaire (on parlera de trend),

La première étape pourrait être de tenter un modèle sans rien.

> 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

La région critique du test est ici (pour un seuil à 5%) l’ensemble des valeurs inférieures à -1.95. Or ici la statistique de test est ici -0.73, on est alors dans la région d’acceptation du test, et on va retenir l’hypothèse de racine unité, i.e. la série n’est pas stationnaire. Mais peut-être faudrait-il juste prendre en compte une constante ?

> 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

Deux statistiques de test sont calculées, ici: la première relative à la racine unité, la seconde relative à la constante.  On observe ici que la statistique de test pour le test de racine unité (-2.3) est ici supérieure à toutes les valeurs critiques associées (données sur la première ligne). On va encore accepter l’hypothèse de racine unité. Mais le modèle était peut-être faux, et peut-être avait-on en fait en tendance linéaire ?

> 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

On obtient cette fois trois statistique, la première relative au test de racine unité, et les deux suivantes aux tests sur la constante, et sur la tendance (la pente de l’ajustement linéaire). Là encore, la valeur de test (-1.98) excède les valeurs critiques: la p-value serait ici supérieure à 10%. On va là encore accepter l’hypothèse de racine unité.
Mais la modélisation autorégressive à l’ordre 1 était peut-être elle aussi fausse. Aussi, il existe un test de Dickey-Fuller augmenté. L’idée est de considérer, de manière beaucoup plus générale, un processus autorégressif à un ordre supérieur,
http://freakonometrics.blog.free.fr/public/perso6/adf-2.gif
Comme auparavant, on peut réécrire ce modèle en fonction des variations (on parlera de modèle à correction d’erreur)
http://freakonometrics.blog.free.fr/public/perso6/acf-29.gif
où on construit les nouveaux coefficients par une relation de récurrence, du genre http://freakonometrics.blog.free.fr/public/perso6/adf-6.gif. Bref, à nouveau, on est ramené à tester http://freakonometrics.blog.free.fr/public/perso6/acf-22.gif ou http://freakonometrics.blog.free.fr/public/perso6/adf-26.gif . Et là encore, on a le choix sur la tendance.
Si on commence par supposer qu’il n’y a pas de tendance,

> df=ur.df(X,type="none",lags=2)
> 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.86738 -0.53887 -0.02009  0.67058  2.45970

Coefficients:
Estimate Std. Error t value Pr(>|t|)
z.lag.1     -0.005168   0.007389  -0.699    0.485
z.diff.lag1 -0.029619   0.065289  -0.454    0.650
z.diff.lag2 -0.037856   0.065436  -0.579    0.563

Residual standard error: 0.9685 on 234 degrees of freedom
Multiple R-squared: 0.004704,	Adjusted R-squared: -0.008056
F-statistic: 0.3687 on 3 and 234 DF,  p-value: 0.7757

Value of test-statistic is: -0.6994

Critical values for test statistics:
1pct  5pct 10pct
tau1 -2.58 -1.95 -1.62

la valeur de la statistique de test excède là encore les valeurs critiques, i.e. la p-value dépasse 10%. Si on rajoute une constante,

> df=ur.df(X,type="drift",lags=2)
> 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.91609 -0.56702  0.01025  0.62109  2.43970

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.31105    0.13332   2.333   0.0205 *
z.lag.1     -0.03744    0.01565  -2.392   0.0175 *
z.diff.lag1 -0.01917    0.06483  -0.296   0.7677
z.diff.lag2 -0.02794    0.06496  -0.430   0.6676
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9594 on 233 degrees of freedom
Multiple R-squared: 0.02664,	Adjusted R-squared: 0.0141
F-statistic: 2.125 on 3 and 233 DF,  p-value: 0.09774

Value of test-statistic is: -2.3924 2.9709

Critical values for test statistics:
1pct  5pct 10pct
tau2 -3.46 -2.88 -2.57
phi1  6.52  4.63  3.81

on valide encore l’hypothèse de racine unité et de même avec une tendance linéaire,

> df=ur.df(X,type="trend",lags=2)
> 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.85835 -0.58826 -0.00867  0.61407  2.47280

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.3530591  0.1522373   2.319   0.0213 *
z.lag.1     -0.0338831  0.0168523  -2.011   0.0455 *
tt          -0.0005674  0.0009879  -0.574   0.5663
z.diff.lag1 -0.0237595  0.0654158  -0.363   0.7168
z.diff.lag2 -0.0328215  0.0656102  -0.500   0.6174
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9608 on 232 degrees of freedom
Multiple R-squared: 0.02802,	Adjusted R-squared: 0.01126
F-statistic: 1.672 on 4 and 232 DF,  p-value: 0.1573

Value of test-statistic is: -2.0106 2.0849 3.0185

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

On notera toutefois que ces trois modèles nous suggèrent de ne pas retenir autant de retard, qui ne semblent pas significatifs.
On peut bien entendu faire ces tests sur de vraies données par exemple le niveau du Nil

> library(datasets)
> NILE=Nile

ou encore les taux d’intérêt américains,

> base=read.table(
+ "http://freakonometrics.free.fr/basedata.txt",
+ header=TRUE)
> Y=base[,"R"]
> Y=Y[(base$yr>=1960)&(base$yr<=1996.25)]
> TAUX=ts(Y,frequency = 4, start = c(1960, 1))

Par exemple, sur cette dernière, on rejette l’hypothèse de stationnarité

>  df=ur.df(y=TAUX,lags=3,type="drift")
>  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
-3.1982 -0.2947 -0.0629  0.3524  3.1899

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.40609    0.16706   2.431 0.016362 *
z.lag.1     -0.06339    0.02500  -2.535 0.012354 *
z.diff.lag1  0.32613    0.08145   4.004 0.000102 ***
z.diff.lag2 -0.31102    0.08027  -3.875 0.000165 ***
z.diff.lag3  0.28712    0.08103   3.543 0.000541 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.7805 on 137 degrees of freedom
Multiple R-squared: 0.2051,	Adjusted R-squared: 0.1819
F-statistic: 8.839 on 4 and 137 DF,  p-value: 2.227e-06

Value of test-statistic is: -2.5354 3.2457

Critical values for test statistics:
1pct  5pct 10pct
tau2 -3.46 -2.88 -2.57
phi1  6.52  4.63  3.81

Mais toutes sortes d’autres tests (plus robustes) peuvent être faits. Les plus connus sont le test de Philips-Perron et le test dit KPSS. Pour ce dernier, il faut spécifier s’il l’on suppose que la série est de moyenne constante, ou si une tendance doit être prise en compte. Si l’on suppose une constante non-nulle, mais pas de tendance, on utilise

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

alors que pour prendre en compte une tendance

> 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

Derrière se cache un teste du multiplicateur de Lagrange. L’hypothèse nulle correspond à l’absence de racine unité: plus la statistique de test est grande, plus on s’éloigne de la stationnarité (hypothèse nulle). Avec ces deux tests, on rejette là encore l’hypothèse de stationnarité de notre marche aléatoire simulée.

Pour le test de Philips-Perron, on a un test de type Dickey-Fuller,

> PP.test(X)

Phillips-Perron Unit Root Test

data:  X
Dickey-Fuller = -2.0116, Trunc lag parameter = 4, p-value = 0.571

Là encore, la p-value nous recommande de valider l’hypothèse de non-stationnarité. A titre de comparaison, si on avait travaillé sur la série différenciée, on aurait accepté l’hypothèse de stationnarité

> PP.test(diff(X))

Phillips-Perron Unit Root Test

data:  diff(X)
Dickey-Fuller = -15.9522, Trunc lag parameter = 4, p-value = 0.01

De l’hebdomadaire au mensuel

Pour le devoir de série temporelles, les données fournies par Google Insight sont hebdomadaires, ce qui peut rendre la modélisation compliquée. Comme on l’a évoqué en fin de cours, il peut être plus simple de travailler sur des données mensuelles. La petite fonction suivante permet de transformer les données pour avoir des données mensuelles (avec des moyennes par mois, ce qui fait qu’un mois de 28 jours et un mois de 31 jours sont comparables). Histoire d’illustrer, prenons un mot au hasard, disons…. épilation. On commence par sauver le fichier csv, et on le lit sous R,

EPILATION=read.table("/home/charpentier/report-epilation.csv", skip=4,header=TRUE,sep=",",nrows=426)

La petite fonction suivant va nous aider à convertir la base en une série temporelle mensuelle,

H2M=function(BASE){ 
X=BASE[,2] 
date1=substr(as.character(BASE$Semaine),1,10) 
date2=substr(as.character(BASE$Semaine),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(Z,start=c(2004,1),frequency=12) 
return(Zts)}

Si on travaille sur la série hebdomadaire, on a la série suivante

hebdo=ts(EPILATION$épilation,start=2004, frequency=52)

La fonction d’autocorrélation est alors

acf(hebdo,160)

Maintenant, on peut travailler sur les données mensuelles. On utilise alors

mensuel=H2M(EPILATION)

Le graphique est alors le suivant

La fonction d’autocorrélation est cette fois la suivante

acf(mensuel,40)

On retrouve le comportement cyclique, avec une saisonnalité annuelle. Mais avec 12 retards, on devrait avoir des modèles plus simples qu’avec 52 retards. Bref, il peut être plus simple de travailler sur des données mensuelles qu’hebdomadaires. Mais chacun est libre de la série qu’il ou elle analysera…

Simulation de séries temporelles

Un billet rapide pour reprendre le code tapé en cours, la semaine passée. Considérons  un processus autorégressif d’ordre 1,  où  est un bruit blanc, stationnaire, i.e.  appartient à l’intervalle . Le code pour simuler un tel processus est

n=1000
bruit=rnorm(n)
phi1= .85
X=rep(NA,n)
X[1]=0
for(t in 2:n){X[t]=phi1*X[t-1]+bruit[t]}
plot(acf(X),lwd=5,col='blue')
plot(pacf(X),lwd=5,col='blue')

ou avec un autocorrélation au premier ordre négative,

phi1= -0.7

On peut aussi regarder un processus autorégressif au second ordre,

sur la figure ci-dessous (avec en haut à gauche le triangle de stationnarité du couple de paramètres).

phi1=  0.3
phi2=  0.5
X=rep(NA,n)
X[1:2]=0
for(t in 3:n){
X[t]=phi1*X[t-1]+phi2*X[t-2]+bruit[t]}

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

theta1=  .8
X=rep(NA,n)
X[1]=0
for(t in 2:n){
X[t]=bruit[t]+theta1*bruit[t-1]}

ou une moyenne mobile du second ordre,

theta1= -.6
theta2=  .5
X=rep(NA,n)
X[1:2]=0
for(t in 3:n){
X[t]=bruit[t]+theta1*bruit[t-1]+
theta2*bruit[t-2]}

 

Inference and autoregressive processes

Consider a (stationary) autoregressive process, say of order 2,

for some white noise  with variance . Here is a code to generate such a process,

> phi1=.5
> phi2=-.4
> sigma=1.5
> set.seed(1)
> n=240
> WN=rnorm(n,sd=sigma)
> Z=rep(NA,n)
> Z[1:2]=rnorm(2,0,1)
> for(t in 3:n){Z[t]=phi1*Z[t-1]+phi2*Z[t-2]+WN[t]}

Here, we have to estimate two sets of parameters: the autoregressive coefficients, and the variance of the innovation process . There are (at least) three techniques to estimate those parameters.

  • using least square regression

A natural idea is to see here a regression model, and thus, if we consider a matrix formulation,

Here we can run (conditional) ordinary least squares estimation,

> base=data.frame(Y=Z[3:n],X1=Z[2:(n-1)],X2=Z[1:(n-2)])
> regression=lm(Y~0+X1+X2,data=base)
> summary(regression)

Call:
lm(formula = Y ~ 0 + X1 + X2, data = base)

Residuals:
Min      1Q  Median      3Q     Max
-4.3491 -0.8890 -0.0762  0.9601  3.6105

Coefficients:
Estimate Std. Error t value Pr(>|t|)
X1  0.45107    0.05924   7.615 6.34e-13 ***
X2 -0.41454    0.05924  -6.998 2.67e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.449 on 236 degrees of freedom
Multiple R-squared: 0.2561,	Adjusted R-squared: 0.2497
F-statistic: 40.61 on 2 and 236 DF,  p-value: 6.949e-16

> regression$coefficients
X1         X2
0.4510703 -0.4145365
> summary(regression)$sigma
[1] 1.449276
  • using Yule-Walker equations

As we’ve seen in class, we can easily get the following equations for the autocovariance functions,

which can also be written

So we just have to solve a simple linear system of equations. Note that if we divide by the variance, those equations can be written in terms of the autocorrelation functions

The code is the following

> rho1=cor(Z[1:(n-1)],Z[2:n])
> rho2=cor(Z[1:(n-2)],Z[3:n])
> A=matrix(c(1,rho1,rho1,1),2,2)
> b=matrix(c(rho1,rho2),2,1)
> (PHI=solve(A,b))
[,1]
[1,]  0.4517579
[2,] -0.4155920

Now, we need to extract the estimated innovation process, from this set of parameters (note that it could be possible to include the variance term in Yule-Walker equations, to get a three dimensional linear equation)

> estWN=base$Y-(PHI[1]*base$X1+PHI[2]*base$X2)
> sd(estWN)
[1] 1.445706

This estimator is probably not the best one (we can take into account that we’ve lost two degrees of freedom), but as a starting point, let us consider this one.

  • using (conditional) likelihood estimators

Finally, we can assume some distribution for the innovation process. Thestandard model is a Gaussian model, i.e.

In that case, the conditional log likelihood (conditional since we set the first two observations here) is

> CondLogLik=function(A,TS){
+ phi1=A[1];  phi2=A[2]
+ sigma=A[3]	; L=0
+ for(t in 3:length(TS)){
+ L=L+dnorm(TS[t],mean=phi1*TS[t-1]+
+ phi2*TS[t-2],sd=sigma,log=TRUE)}
+ return(-L)}

Now, we can run standard optimization procedures,

> LogL=function(A) CondLogLik(A,TS=Z)
> optim(c(0,0,1),LogL)
$par
[1]  0.4509685 -0.4144938  1.4430930

$value
[1] 425.0164

$counts
function gradient
88       NA

$convergence
[1] 0

$message
NULL

Here, our three estimators are rather close. Actually, if we generate 1,000 time series (of size 240), those are the Box-plots of our three estimators, for the first order autoregressive coefficient

for the second one,

and finally for the standard deviation of the innovation process

All those estimators behave nicely, and are rather close. Note that they all might be biased, but they are consistent (see Davidson and MacKinnon for instance, in their book, for more details).