• Les supports pédagogiques sont déposés au dépôt de (https://github.com/Hamrita/TSFin.git).

• Logiciel statistique: (https://www.r-project.org/)

• IDE: RStudio (https://www.rstudio.com/products/rstudio/download/)

• Un processus \(x_t\) est une série temporelle discrète si \(x_t\) est une variable aléatoire et l’indice \(t\) est dénombrable.

• Une série observée est une réalisation de ce processus stochastique.

• Un processus \(x_t\) est strictement stationnaire si sa distribution est invariante dans le temps. Mathématiquement parlant, \(x_t\) est strictement stationnaire si pour tout indice de temps arbitraire \(\{t_1, t_2,\ldots, t_m\}\), où \(m > 0\), et pour un entier \(k\) fixé, \(F(x_{t_1},\ldots,x_{t_m})=F({x_{t_{1}+k}, \ldots, x_{t_{m}+k}})\).

• Une série temporelle est faiblement stationnaire si les deux premiers moments de \(x_t\) existent et invariants dans le temps. Statistiquement parlant, \(\mathbb{E}(x_t)=m\) et \(Cov(x_t,x_{t+k})=\gamma(k)\), où \(\mathbb{E}()\) est l’espérance mathématique, \(Cov\) est la covariance et \(\gamma(k)\), dite fonction d’auto-covariance d’ordre \(k\) vérifiant \(\gamma(-k)=\gamma(k)\), \(\forall k \in \mathbb{Z}\).

• Une séquence de variables aléatoires indépendantes et identiquement distribuées est strictement stationnaire.

Un processus \(x_t\) est dit linéaire s’il peut s’écrire sous la forme: (représentation de Wald) \[ x_t=m+\sum_{k=-\infty}^{+\infty} \psi_k e_{t-k} \] où \(m\) et \(\psi_k\) sont deux réels avec \(\psi_0=1\) et \(\displaystyle \sum_{k=-\infty}^{+\infty} |\psi_k| < \infty\) et \(\{e_t\}\) est une séquence aléatoire \(iid\) de moyenne nulle et admettant une distribution. Dans la pratique, on s’intéresse aux séries temporelles unilatérales

\[ x_t=m+\sum_{k=0}^{+\infty} \psi_k e_{t-k} \]

La série temporelle linéaire dans l’équation précédente est faiblement stationnaire si nous supposons en outre que \(\mathbb{V}(e_t)=\sigma^2_e\).

Dans ce cas, on a \(\mathbb{E}(x_t)=m\), \(\mathbb{V}(x_t)=\displaystyle \sigma^2_e\sum_{k=0}^\infty \psi_k^2\) et \(\gamma(k)=\sigma^2_e\displaystyle\sum_{i=0}^\infty\psi_i\psi_{i+k}\).

Le modèle ARIMA est un modèle célèbre qui admet cette écriture. Dans ce cas, les coefficients \(\psi_i\) se calculent comme suit:

\[ \psi_i=\theta_i+\sum_{0 \leq k \leq i}\phi_k \psi_{i-k}\;;\qquad 0 \leq i < \max(p,q+1) \] \[ \psi_i=\sum_{0 \leq k \leq i}\phi_k \psi_{i-k}\;;\qquad i \geq \max(p,q+1) \]

Soit \(x_t \sim AR(1)\) tel que \(x_t=0.4 x_{t-1}+e_t\); \(\quad e_t\stackrel{iid}{\sim}BB(0,1)\)

1. Déterminer la représentation de Wald et donner les valeurs de \(\psi_i\) pour \(i=0,1,2,\ldots,10\).
2. En déduire l’expression de \(\gamma(k)\) en fonction de \(k\). Donner les valeurs de \(\gamma(k)\) pour \(k=0,1,\ldots, 10\).

psi_i=ARMAtoMA(ar=0.4,lag.max = 10)  # représentation de Wald
psi_i

 [1] 0.4000000000 0.1600000000 0.0640000000 0.0256000000 0.0102400000
 [6] 0.0040960000 0.0016384000 0.0006553600 0.0002621440 0.0001048576

ARMAacv=function(ar=0,ma=0,n=10){
  p=c(1,ARMAtoMA(ar,ma,n+10000))
  gg=NULL
  for(ii in 1:n) gg[ii]=sum(p[1:10000]*p[(ii+1):(10000+ii)])
  c(sum(p^2),gg)
}
ARMAacv(ar=0.4,n=10)

 [1] 1.1904761905 0.4761904762 0.1904761905 0.0761904762 0.0304761905
 [6] 0.0121904762 0.0048761905 0.0019504762 0.0007801905 0.0003120762
[11] 0.0001248305

# 0.4^(0:10)*25/21

• Tout processus qui ne peut pas s’écrire sous la représentation de Wald est un processus non linéaire.

• Nous commençons par un exemple réel qui présente clairement des caractéristiques non linéaires et nécessite une modélisation non linéaire.

• On considère le cours journalier de Bitcoin contre le dollar américain (BTC/USD) allant de 17-09-2014 à 20-10-2021. On donne l’évolution du cours ajusté (figure à gauche) ainsi que l’évolution du rendement qui est défini par \(r_t=\ln\left(\dfrac{p_t}{p_{t-1}}\right)\) où \(p_t\) est le cours ajusté.

library(quantmod); library(zoo)
btc=getSymbols("BTC-USD", src="yahoo", from="2014-09-17",to="2021-10-20",auto.assign = FALSE)
closedAdj=zoo(btc[,6]); rt=diff(log(closedAdj))
par(mfrow=c(1,2), mar=c(2.5,3.8,1,1))
plot(closedAdj, col=4, xlab=""); plot(rt,col=4, xlab="")

• La figure de l’évolution des prix (à gauche) montre que la moyenne ainsi que la variance évoluent dans le temps tous les deux, donc la série est non stationnaire.

• En transformant la série, \(r_t=\ln(p_t)-\ln(p_{t-1})\), on voit bien que la moyenne semble être stable dans le temps. Tandis que la variance semble être variable dans le temps et de manière asymétrique. Ce qui suggère à son tour que les rendements du Bitcoin ne sont pas linéaires.

par(mfrow=c(2,2), mar=c(2.5,3.8,1,1))
acf(closedAdj, 40,na.action = na.pass, col=4, xlab=""); acf(rt, 40,na.action = na.pass, col=4, xlab="")
pacf(closedAdj,40,na.action = na.pass,col=4, xlab=""); pacf(rt,40,na.action = na.pass,col=4, xlab="")

• grande variété des séries utilisées (prix d’action, taux d’intérêt, taux de change etc.), importance de la fréquence d’observation (seconde, minute, heure, jour, etc), disponibilité d’échantillons de très grande taille.

• existence de régularités statistiques (“faits stylisés”) communes à un très grand nombre de séries financières et difficiles à reproduire artificiellement à partir de modèles stochastiques. Mandelbrot (1963).

  • Non stationnarité des prix \(p_t\)
  • Possible stationnarité des rendements.
  • Regroupement des extrêmes (volatility clustering): On observe empiriquement que de fortes variations des rendements sont généralement suivies de fortes variations. Ce type de phénomène remet en cause l’hypothèse d’homoscédasticité.
  • Non corrélation des rendements mais auto-corrélation des carrés.
  • Queues de distribution épaisses: la distribution des résidus demeure leptokurtique.
  • Asymétrie: Les baisses du cours génèrent plus de volatilité que les hausses de même amplitude (effet de levier).
  • Saisonnalité: Les rendements présentent de nombreux phénomènes de saisonnalité (effets week end, effet janvier, etc..)

• La flexibilité des modèles non linéaires dans l’ajustement des données peut rencontrer le problème de trouver une structure fallacieuse (spurious) dans une série chronologique donnée.

• Il est donc important de vérifier la nécessité d’utiliser des modèles non linéaires. À cette fin, nous introduisons des tests de linéarité pour les données de séries temporelles. Les tests paramétriques et non paramétriques sont pris en compte.

• Tests paramétriques:

  • Le test de multiplicateur de Lagrange ( Engle (1982) )
  • Le test RESET (Regression error specification test) ( Ramsey (1969), Keenan (1985) )

• Tests non paramétriques

  • Le test McLeod-Li ( McLeod and Li (1983) )
  • Le test BDS ( Broock et al. (1996) )
  • Test de Peña-Rodríguez (Peña and Rodríguez 2002, 2006)

Tests paramétriques

Le test LM: Il s’agit de tester l’effet ARCH (héteroscédasticité conditionnelle) proposé par Engle (1982). L’hypothèse nulle est: \(H_0: \alpha_1 = \alpha_2 = \ldots = \alpha_m = 0\) du modèle: \(e_t^2=\alpha_0+\alpha_1 e_{t-1}^2 + \alpha_2 e_{t-2}^2 + \ldots = \alpha_m e_{t-m}^2 + a_t\), \(t=1,2,\ldots, T\).

\(e_t\) est la série des résidus obtenue suite à un modèle linéaire (ARIMA, régression linéaire).

\[LM=T R^2 \stackrel{H_0}{\sim}\chi^2(m) \text{ où }R^2 \text{ est le coefficient de détermination.} \]

Si \(LM < \chi^2_{\alpha}(m)\), l’hypothèse nulle est acceptée; absence d’effet ARCH.

Remarques:

• Ce test est équivalent au test de Fisher; \[ F=\dfrac{(SCR_0-SCR_1)/m}{SCR_1/(T-2m-1)}\stackrel{H_0}{\sim}F(m,T-2m-1) \] où \(SCR_0\) et \(SCR_1\) sont les sommes des carrés résiduelles sous \(H_0\) et \(H_1\) respectivement.

• On accepte l’hypothèse nulle si p.value est supérieure à \(\alpha \%\).

archLM=function(x,p=1, disp=T){
  # x: time series;  p: lag order; disp: display or not the result
  mat=embed(x^2,(p+1))
  reg=summary(lm(mat[,1]~mat[,-1]))
  r2=reg$r.squared
  statistic=r2*length(resid(reg))
  cv=qchisq(0.95,p)
  pvalue=1-pchisq(statistic,p)
  if(disp){
  cat("================================================\n")
  cat("=======           ARCH LM test         =========\n")
  cat("================================================\n")
  cat("====   H0: no ARCH effects             =========\n\n")
  
  cat(" Chi-squared: ", statistic, "  df: ", p,  "  p.value: ", pvalue, 
      "  critical value 5%: ", cv, " \n")
  cat(" F-statistic", reg$fstatistic[1]," df:  ", reg$fstatistic[-1],"p.value: ",
      1-pf(reg$fstatistic[1], reg$fstatistic[2],reg$fstatistic[3]), 
      " critical value 5%: ", qf(0.95,reg$fstatistic[2], reg$fstatistic[3]), "\n")
  }
  return(invisible(list("Chi-squared"=statistic, "p-value"=pvalue,
              "df"=p)))
}

Appliquons le test LM sur les rendements du Bitcoin (en supposant que \(r_t=m+\epsilon_t\))

et=rt-mean(rt,na.rm=T)
archLM(et, 12)

================================================
=======           ARCH LM test         =========
================================================
====   H0: no ARCH effects             =========

 Chi-squared:  79.22112   df:  12   p.value:  5.809353e-12   critical value 5%:  21.02607  
 F-statistic 6.778986  df:   12 2533 p.value:  3.736567e-12  critical value 5%:  1.755988 

Le résultat montre bien l’existence d’un effet ARCH sur les rendements du Bitcoin (p.value presque nulle).

RESET Test

Ramsey (1969) a proposé un test de spécification pour les modèles de régression linéaire. Ce test est aussi applicable aux modèles \(AR\). Les étapes du test sont décrites comme suit:

• Estimation, par MCO, le modèle: \(x_t=\phi_0+\phi_1 x_{t-1}+\ldots + \phi_p x_{t-p}+\varepsilon_t\), on obtient, alors \(\mathbf{\widehat{\phi}}=(\widehat{\phi}_0,\widehat{\phi}_1,\ldots, \widehat{\phi}_p)\), \(\widehat{x}_t\), \(e_t=\widehat{\varepsilon}_t=x_t-\widehat{x}_t\) et \(SCR_0=\sum e_t^2\).

• Estimation, par MCO, le modèle: \(x_t=\alpha_0+\alpha_1 x_{t-1}+\ldots + \alpha_p x_{t-p}+ \beta_{\color{red}{1}} \widehat{x}_t^{\color{red}{2}}+\beta_{\color{red}{2}} \widehat{x}_t^{\color{red}{3}}+ \ldots +\beta_{\color{red}{s}} \widehat{x}_t^{\color{red}{s+1}}+\nu_t\), pour \(s \geq 1\). En déduire \(SCR_1=\sum \widehat{\nu}_t^2\). Le modèle est non linéaire, si \(H_0: \alpha_0=\alpha_1=\ldots=\beta_1=\beta_s=0\), donc on peut faire recours au test du Fisher.

• Calcul du statistique \(F=\dfrac{(SCR_0-SCR_1)/(s+p+1)}{SCR_1/(T-2p-s-1)}\stackrel{H_0}{\sim}F(s+p+1,T-2p-s-1)\)

• Remarque: Parce que les variables \(\widehat{x}_j\) pour \(j=2, \ldots, s+1\) ont tendance à être fortement corrolées avec \((x_{t-1},x_{t-2}, \ldots, x_{t-p})\) et entre elles, les composantes principales de \((\widehat{x}_t^2, \ldots, \widehat{x}_t^{s+1})\) qui ne sont pas colinéaires avec \((x_{t-1},x_{t-2}, \ldots, x_{t-p})\) sont souvent utilisées dans l’ajustement de la deuxième équation.

reset.test=function(x,p=1,s=1){
  # Estimate of AR(p)
  xx=embed(x,(p+1))
  mod1=lm(xx[,1]~xx[,-1])
  scr0=sum(mod1$residuals^2)
  xx.hat=mod1$fitted.values
  et=mod1$residuals;  TT=length(xx.hat)
  # Estimate of model 2
  xxx=matrix(0,nr=TT, nc=s)
  for(ii in 1:s)  xxx[,ii]=xx.hat^(s+1)
  X=cbind(xx[,-1],xxx)
  mod2=lm(et~X);  scr1=sum(mod2$residuals^2)
  df2=TT-p-s-1 ; df1=s+p+1
  F_stat=df2/df1*(scr0-scr1)/scr1
  pval=pf(F_stat,df1,df2,lower.tail = F)
  cval=qf(0.95,df1,df2)
  dec=ifelse(pval < 0.05, "H0 is rejected", "H0 is accepted")
  cat("=======================================\n")
  cat("        RESET test                     \n")
  cat("=======================================\n\n")
  cat("      H0: the model is linear          \n\n")
  cat("F-statistic: ", F_stat, "df: ", c(df1,df2), " p-value: ", pval, "CV 5%: ", cval,"\n")
  cat("Decision: ", dec , "\n")
  
  return(invisible(list(Fstatistic=F_stat,df=c(df1,df2),p.value=pval)))
}

reset.test(na.omit(rt), p=2,1)

=======================================
        RESET test                     
=======================================

      H0: the model is linear          

F-statistic:  1.298354 df:  4 2577  p-value:  0.2683325 CV 5%:  2.375381 
Decision:  H0 is accepted 

Keenan Test

• Keenan (1985) a proposé un test de non-linéarité pour les séries temporelles qui utilise uniquement \(\widehat{x}^2_t\) et modifie la deuxième étape du test RESET pour éviter la multi-colinéarité entre \(\widehat{x}^2_t\) et \((x_{t-1},x_{t-2}, \ldots, x_{t-p})\).
• Plus précisément, la deuxième régression linéaire est divisée en deux étapes. En premier lieu, on supprime la dépendance linéaire de \(\widehat{x}^2_t\) sur \((x_{t-1},x_{t-2}, \ldots, x_{t-p})\) en ajustant la régression \(\widehat{x}^2_t=\alpha_0+\alpha_1 x_{t-1}+\ldots+\alpha_p x_{t-p}+v_t\) et déduire les résidus \(\widehat{v}_t\).
• En second lieu, on considère le modèle \(e_t = \gamma \widehat{v}_t + u_t\) pour obtenir la somme des carrées \(SCR_0=\sum\widehat{u}^2_t\) afin de tester l’hypothèse nulle \(H0: \gamma =0\)

keenan.test=function(x,p=1){
  # Estimate of AR(p)
  xx=embed(x,(p+1))
  mod1=lm(xx[,1]~xx[,-1])
  xx.hat=mod1$fitted.values
  et=mod1$residuals
  # Estimate of model 2
  # (i)
  mod21=lm(xx.hat^2~xx[,-1])
  vt=mod21$residuals
  mod22=lm(et ~ -1+vt)
  F_stat=summary(mod22)$fstatistic[1]
  df=summary(mod22)$fstatistic[-1]
  pval=pf(F_stat,df[1],df[2],lower.tail = F)
  cval=qf(0.95,df[1],df[2])
  dec=ifelse(pval < 0.05, "H0 is rejected", "H0 is accepted")
  cat("=======================================\n")
  cat("        Keenan test                    \n")
  cat("=======================================\n")
  cat("      H0: the model is linear          \n\n")
  print(summary(mod22)$coefficients)
  cat("\n\nF-statistic: ", F_stat, "df: ", c(df), " p-value: ", pval, "CV 5%: ", cval,"\n")
  cat("Decision: ", dec , "\n")
  
  return(invisible(list(Fstatistic=F_stat,df=c(df),p.value=pval)))
}

keenan.test(na.omit(rt),2)

=======================================
        Keenan test                    
=======================================
      H0: the model is linear          

   Estimate Std. Error  t value   Pr(>|t|)
vt  565.753    248.112 2.280232 0.02267514


F-statistic:  5.19946 df:  1 2580  p-value:  0.02267514 CV 5%:  3.845066 
Decision:  H0 is rejected 

Le test de McLeod-Li

• Une statistique de type test-portemanteau, basée sur la fonction d’auto-corrélation des carrés résiduels obtenus à partir d’un modèle ARMA, a été proposée par McLeod and Li (1983). L’idée est d’appliquer la statistique de Ljung-Box aux carrés des résidus d’un modèle \(ARMA(p,q)\) pour vérifier l’inadéquation du modèle. Par conséquent, la statistique de test est \[ Q(m)=n(n+2)\sum_{i=1}^n\dfrac{\widehat{\rho}^2_i(e^2_i)}{n-i} \] où \(n\) est la taille de l’échantillon, \(m\) est un nombre correctement choisi d’auto-corrélations utilisées dans le test et \(e_i\) les résidus du modèle ARMA.
• Sous \(H_0\), \(Q(m) \sim \chi^2(m-p-q)\).

arma_pq=arima(rt,c(1,0,1))
resid2=arma_pq$residuals^2
m=1L:7L; Q=NULL ; pval=NULL
for (i in m) { Q[i]=Box.test(resid2,i,"Ljung-Box")$statistic
pval[i]= Box.test(resid2,i,"Ljung-Box")$p.value
}
tab=rbind(m,Q,p.value=pval)
print(tab, digits=4)

             [,1]      [,2]      [,3]      [,4]      [,5]      [,6]     [,7]
m       1.000e+00 2.000e+00 3.000e+00 4.000e+00 5.000e+00 6.000e+00 7.00e+00
Q       3.477e+01 3.973e+01 4.203e+01 5.550e+01 6.018e+01 6.433e+01 8.83e+01
p.value 3.719e-09 2.358e-09 3.954e-09 2.554e-11 1.116e-11 5.898e-12 2.22e-16

• Le test BDS : Le test BDS (Broock et al. (1996)), développé dans le cadre de la théorie du chaos, est l’un des tests de non-linéarité les plus populaires.
• La statistique BDS est basée sur l’intégrale des corrélations. Etant donné une série temporelle \(\{x_t\}\), et en posant \(x_t^m=\{x_t, x_{t-1},\ldots,x_{t-m+1}\}\), l’intégrale des corrélations est donnée par: \[ C_{m,T}(\epsilon)=\sum_{t < s}I_{\epsilon}x^m_t x^s_t \left\{\dfrac{2}{T_m(T_m-1)} \right\} \] où \(T_m=T−(m−1)\) et \(I_{\epsilon}x^m_tx^m_s\) est une fonction indicatrice égale à un si \(||x^m_t-x^s_t|| < \epsilon\) et égale à zéro sinon avec \(||.||\) désigne la norme supérieure.
• La statistique de BDS est donnée par: \[ W_{m \epsilon}=\sqrt{T}\dfrac{C_{m,T}(\epsilon)-C_{1,T}(\epsilon)^m}{s_{m,T}} \] où \(s_{m,T}\) est l’écart-type estimé.
• Sous l’hypothèse nulle, pour laquelle la série \(\{x_t\}\) est considérée \(iid\), Broock et al. (1996) ont montré que \(W_{m \epsilon}\) converge vers la loi normale standard.
• Dans la pratique, le test BDS est appliqué aux résidus suite à une estimation d’un modèle linéaire.

mod=arima(rt,c(1,0,1))
residd=c(mod$residuals) ; 
fNonlinear::bdsTest(na.omit(residd))@test

$data.name
[1] "na.omit(residd)"

$statistic
eps[1] m=2 eps[1] m=3 eps[2] m=2 eps[2] m=3 eps[3] m=2 eps[3] m=3 eps[4] m=2 
 10.947219  14.762387   9.236784  12.056526   8.272856  10.621231   7.758839 
eps[4] m=3 
  9.459227 

$p.value
  eps[1] m=2   eps[1] m=3   eps[2] m=2   eps[2] m=3   eps[3] m=2   eps[3] m=3 
6.851911e-28 2.560545e-49 2.540183e-20 1.791830e-33 1.307873e-16 2.374085e-26 
  eps[4] m=2   eps[4] m=3 
8.571025e-15 3.102280e-21 

$parameter
Max Embedding Dimension                  eps[1]                  eps[2] 
             3.00000000              0.01974641              0.03949281 
                 eps[3]                  eps[4] 
             0.05923922              0.07898562 

Test de Peña-Rodríguez (Peña and Rodríguez 2002, 2006)

• Peña et Rodríguez (2002) ont proposé une statistique de test de Portemanteau qui peut être utilisée pour la vérification d’un modèle linéairement ajusté, y compris la non-linéarité dans les résidus.
• Soit \(e_t\) les résidus suite à un modèle linéaire. Soit \(z_t\) une fonction de \(e_t\); \(z_t=e^2_t\) ou \(z_t=|e_t|\). L’auto-corrélation empirique de retard \(k\) de \(z_t\) est: \[ \widehat{\rho}_k=\dfrac{\displaystyle{\sum_{i=k+1}^T(z_{t-i}}-\overline{z})(z_t-\overline{z})}{\displaystyle{\sum_{i=1}^T}(z_t-\overline{z})^2} \]
• Pour un entier positif donné, \(m\), l’hypothèse nulle du test de Peña-Rodríguez est \(H_0: \rho=\rho_2=\ldots=\rho_m\). La statistique est \(\widehat{D}_m=T\left(1-\left|\widehat{R}_m \right|^{1/m} \right)\) où \[ \widehat{R}_m=\left(\begin{array}{cccc} 1 & \widehat{\rho}_1 & \ldots & \widehat{\rho}_m\\ \widehat{\rho}_1 & 1 & \ldots & \widehat{\rho}_{m-1}\\ \cdots & \cdots & \ddots & \vdots\\ \widehat{\rho}_m & \widehat{\rho}_{m-1} & \ldots & 1 \end{array}\right) \]

• En utilisant l’idée de pseudo-vraisemblance, (Peña and Rodríguez 2006) ont modifié la statistique du test: \[ \widehat{D}^*_m=-\dfrac{T}{m+1}\log\left(|\widehat{R}| \right) \] qui est distribuée asymptotiquement comme un mélange de \(m\) variables aléatoires indépendantes \(\chi^2(1)\). Cependant, les auteurs ont dérivé deux approximations pour simplifier le calcul. Ils ont aboutit au résultat \(N\widehat{D}^*_m \sim N(0,1)\) où \(N\widehat{D}^*_m\) est une approximation à \(\widehat{D}^*_m\) (Eq 10 de Peña and Rodríguez (2006)).

• Sous , on peut faire recours à la fonction PRnd du package NTS.

NTS::PRnd(na.omit(residd))

lambda:  3.972973 
ND-stat & p-value  0.1135843 0.9095673 

Références