{"word_count": 530, "character_count": 4311}
Yesterday, we've seen how inference for  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  time seriesfor some white noise .&#62; theta1=.25&#62; theta2=.7&#62; n=1000&#62; set.seed(1)&#62; e=rnorm(n)&#62; Z=rep(0,n)&#62; for(t in 3:n) Z[t]=e[t]+theta1*e[t-1]+theta2*e[t-2]&#62; Z=Z[800:1000]&#62; plot(Z,type="l")Using the empirical autocorrelationsThe first idea might be to use the first two (empirical) autocorrelations (the two that are supposed to be - theoretically - non null).with  when . We also have the following relationship on the variance of the processWith those three equations, for three unknown parameters, ,  and , we simply have to solve (numerically) that system of equations,&#62; v=c(as.numeric(acf(Z)$acf[2:3]),var(Z))&#62; v[1] 0.1658760 0.3823053 1.6379498&#62; library(rootSolve)&#62; 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))}&#62; 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-09We 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 techniquesWe can try something quite different here. The problem we have is that we do not observe the noise , we only observe our series . But we can try to rebuild that series (call it  since we're not sure it will be a reconstruction of the noise). As suggested in Box &#38; Jenkins (1967), assume that the first two values are null. And then, useand then, we can use least square techniquesThe code will be&#62; 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&#62; optim(par=c(.1,.1),V)$par[1] 0.2751667 0.6723909$value[1] 225.8104$countsfunction gradient       77       NA $convergence[1] 0$messageNULLwhich 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,&#62; plot(e[800:1000],col="blue",type="l")&#62; theta1=0.2751667&#62; theta2=0.6723909&#62; u=rep(0,length(Z))&#62; for(t in 3:length(Z)) u[t]=Z[t]-theta1*u[t-1]-theta2*u[t-2]&#62; lines(1:201,u,col="red")So far, so good. And if we look at 1,000 samples, we getIt 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 techniqueAnd 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  is Gaussian, with a simple variance matrix (since a lot of elements will be null),&#62; library(mnormt)&#62; 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))}&#62; LogL=function(A) -GlobalLogLik(A,TS=Z)&#62; optim(c(.1,.1,1),LogL)$par[1] 0.2584144 0.6826530 1.0669820$value[1] 298.8699$countsfunction gradient       86       NA $convergence[1] 0$messageNULLHere, 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  time series.