More on trends with autocorrelation

In the previous autocorrelation post I described a general multiple regression method with AR(n) noise models. It showed how jointly regressing a trend variable with lagged data led tp a polynomial equation in the coefficients which could be rapidly solved by a Newton-Raphson iteration.

Hu McCulloch noted a similarity with Cochrane-Orcutt. For AR(1) it is similar, though C-O descriptions are usually vague on how iterations should be done.

In this post I want to spell out the derivative calculation in the iteration. This showed me how to write the code more systematically, so I'll include below a general AR(n) version.

In the earlier post I got to the stage where
E_kimja_ia_jb_mb_k is to be minimised, where a are the AR(n) coefficients, b the (2) linear regression coefficients, and
E = E_kimj=Lⁱ_pq yX_qk L^j_pr yX_rm
See the earlier post for further definitions and the use of subscripts and the summation convention. For consistency with the computer program, I've changed the order of subscripts on E, which has symmetry between k and m, and also i and j.
Minimisation is done by differentiating wrt a and b, remembering that the first coefficient is fixed at -1. So
E_kimja_ia_jb_m =0
E_kimja_jb_kb_m =0
I've shown free (unpaired) suffices in red. In Newton Raphson terms, this is the equation f(c)=0, where c=(b,a). And we have to iteratively solve for
f(c)+f'(c)δc=0. Writing this out:
E_kimja_ia_jb_m + E_kimja_ia_j δb_m + 2*E_kimja_ib_m δa_j = 0
E_kimja_jb_kb_m + 2*E_kimja_jb_k δb_m + E_kimjb_kb_m δa_j = 0
So here is the program to implement this:

# A program written by Nick Stokes to get trend of a time series allowing for AR(n) autocorrelation of noise

k=5;  # Setting p for AR(p);  k=p+1
# define a function for multi-index array multiplication
mul=function(x,y,n){apply(x,(1:length(dim(x)))[-n],"%*%",y)}

#  Read Tamino's file, which I also use for HADCRUT
v=read.csv("AllFit.csv")
#  Extract the wanted time series
g=window(ts(v[,8],s=1950,f=12),start=c(1994,12),end=2010)
n=length(g); a=NULL

# Make X (yX), the regression matrix with y added. Time in century units
X=cbind(g,1,(1:n-n/2-.5)/1200)
# Add lagged X

i=k:n; a=Xi=NULL
#  Make matrix of regressors with lags
for(j in 1:k-1)Xi=cbind(Xi,X[i-j,])
# Form all products
E=t(Xi)%*%Xi
#  Note - that's the end of time series. Now just small arrays
#  Now get the OLS regression coefs as a starting point for Newton
df=solve(E[2:3,2:3])
#  We're fitting the AR coefs - they start as 0
b = c( df %*% E[2:3,1],rep(0,k-1))
# Set up E as 4-index array
dim(E)=c(3,k,3,k)
#  Now the Newton Raphson iteration which fits the coefficients
for(i in 1:6){
#  b is the vector of unknowns. Split into a0,b0
  b0=c(-1,b[1:2])
  a0=c(-1,b[2:k+1])
# Make the second derivative blocks  
  c2=mul(mul(E,b0,3),b0,1)
  c1=mul(mul(E,a0,2),2*b0,1)
  c0=mul(mul(E,a0,4),a0,2)
  j = -c(1,4)  # For removing rows/columns corresp to fixed coefs (-1)
#  Make f and f' for Newton step
  f=rbind(c0%*%b0,c2%*%a0)
  df=rbind(cbind(c0,c1),cbind(t(c1),c2))[j,j]
#  s is the sum of squares being minimised
  s=sum(f*c(b0,a0))/2 
# Store items for printing
  a=rbind(a,c(b0[3],b0[3]/sqrt(s*df[2,2]/(n-4)),s))
#  Invert the derivative and increment b
  df=solve(df)
  db=-df%*%f[j]
  b=b+db
}
colnames(a)=c("Slope","  t-statistic"," Sum Squares")
print(a)

Synthetic tests

I had promised in the last post to check with series from Foster and Rahmstorf, where they do an enhanced Ar(1) check. But I see that they have added some harmonics of the annual to remove some residual periodicity. I didn't want to deal with that. So I tried generating some synthetic AR(1) data.
I encountered some difficulty using the R routine arima.sim, in that the rho values did not quite correspond to the AR() parameters. The effect seemed to depend on the length of sequence, approaching the expected value for long series, so it may be an end effect. However, I found that they did correspond to ρ directly estimated. So I can use the series.

So here is the first test. I generated 1000 series of M=100 terms using noise=arima.sim(list(ar=c(.8)),M) ts=(im/1200 + gg) The directly calculated ρ turned out to be 0.787&plusmin;0.063.

....more to come