# Inference for MA(q) Time Series

Yesterday, we’ve seen how inference for $AR(p)$ time series was possible.  I started  with that one because it is actually the simple case. For instance, we can use ordinary least squares. There might be some possible bias (see e.g. White (1961)), but asymptotically, estimators are fine (consistent, with asymptotic normality). But when the noise is (auto)correlated, then it is more complex. So, consider here some $MA(2)$ time series

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

for some white noise $(\varepsilon_t)$.

> theta1=.25
> theta2=.7
> n=1000
> set.seed(1)
> e=rnorm(n)
> Z=rep(0,n)
> for(t in 3:n) Z[t]=e[t]+theta1*e[t-1]+theta2*e[t-2]
> Z=Z[800:1000]
> plot(Z,type="l")

• Using the empirical autocorrelations

The first idea might be to use the first two (empirical) autocorrelations (the two that are supposed to be – theoretically – non null).

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

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

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

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

With those three equations, for three unknown parameters, $\theta_1$$\theta_2$ and $\sigma$, we simply have to solve (numerically) that system of equations,

> v=c(as.numeric(acf(Z)$acf[2:3]),var(Z)) > v [1] 0.1658760 0.3823053 1.6379498 > library(rootSolve) > seteq=function(x){ + F1=v[1]-(x[1]+x[1]*x[2])/(1+x[1]^2+x[2]^2) + F2=v[2]-(x[2])/(1+x[1]^2+x[2]^2) + F3=v[3]-(1+x[1]^2+x[2]^2)*x[3]^2 + return(c(F1,F2,F3))} > multiroot(f=seteq,start=c(.1,.1,1))$root
[1] 0.1400579 0.4766699 1.1461636

$f.root [1] 7.876355e-10 4.188458e-09 -2.839977e-09$iter
[1] 5

$estim.precis [1] 2.605357e-09 We are a bit far away from the true values, used to generate our sample. And if we consider 1,000 sample (instead of only one), we still have the bias, and a large variance for our three estimators, • Using least square techniques We can try something quite different here. The problem we have is that we do not observe the noise $(\varepsilon_t)$, we only observe our series $(X_t)$. But we can try to rebuild that series (call it $(u_t)$ since we’re not sure it will be a reconstruction of the noise). As suggested in Box & Jenkins (1967), assume that the first two values are null. And then, use $u_t =X_t-\left(\theta_1u_{t-1}+\theta_2un_{t-2}\right)$ and then, we can use least square techniques $\min_{(\theta_1,\theta_2)}\left\{\sum_{i=1}^tu_i^2\right\}$ The code will be > V=function(p){ + theta1=p[1] + theta2=p[2] + u=rep(0,length(Z)) + for(t in 3:length(Z)) u[t]=Z[t]-theta1*u[t-1]-theta2*u[t-2] + return(sum(u^2)) + } If we try to minimize the sum of the squares of the residuals, we get > optim(par=c(.1,.1),V)$par
[1] 0.2751667 0.6723909

$value [1] 225.8104$counts
77       NA

$convergence [1] 0$message
NULL

which is close to the true value. Another good thing is that, if we compare that rebuilt noise with the true one (since we actually have it), then we have the same vector,

> plot(e[800:1000],col="blue",type="l")
> theta1=0.2751667
> theta2=0.6723909
> u=rep(0,length(Z))
> for(t in 3:length(Z)) u[t]=Z[t]-theta1*u[t-1]-theta2*u[t-2]
> lines(1:201,u,col="red")

So far, so good. And if we look at 1,000 samples, we get

It looks like we have some bias here. And since the two estimators should be negatively correlated, one over-estimates, while the other one under-estimates.

• Using the (global) maximum likelihood technique

And a final method might be to use the maximum likelihood technique (globally). Again, if we assume that we have a Gaussian i.i.d noise, then the vector $\boldsymbol{Y}=(Y_1,\cdots,Y_t)$ is Gaussian, with a simple variance matrix (since a lot of elements will be null),

> library(mnormt)
> GlobalLogLik=function(A,TS){
+ n=length(TS)
+ theta1=A[1];  theta2=A[2]
+ sigma=A[3]
+ SIG=matrix(0,n,n)
+ rho=rep(0,n)
+ rho[1]=1
+ rho[2]=(theta1+theta1*theta2)/(1+theta1^2+theta2^2)
+ rho[3]=(theta2)/(1+theta1^2+theta2^2)
+ for(i in 1:n){for(j in 1:n){
+ SIG[i,j]=rho[abs(i-j)+1]}}
+ gamma0=(1+theta1^2+theta2^2)*sigma^2
+ SIG=gamma0*SIG
+ return(dmnorm(TS,rep(0,n),SIG,log=TRUE))}
> LogL=function(A) -GlobalLogLik(A,TS=Z)
> optim(c(.1,.1,1),LogL)
$par [1] 0.2584144 0.6826530 1.0669820$value
[1] 298.8699

$counts function gradient 86 NA$convergence
[1] 0

\$message
NULL

Here, the values that minimize the likelihood are rather close to the ones used to generate our sample. And if we run this algorithm on 1,000 samples, we can see that those estimates are fine,

I could not find other ideas, to estimate those parameters. I guess we can use the partial autocorrelation function, since we have relationships that can be related to Yule-Walker equations for $AR(p)$ time series.

# 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]}