# Why does accounting for autocorrelated residuals barely help parameter estimation in distributed lag models

This problem has been plaguing me for a long time. Basically, I have a distributed lag model $$y_t=\sum_{i=0}^{p} \beta_i x_{t-i} + u_t.$$

The regression problem is a bit misspecified, so I end up with autocorrelated errors $$u_t=\alpha u_{t-1}+\epsilon_t.$$

Due to the high levels of autocorrelation in my model, the estimation of $$\beta$$ should be very inefficient, but when I correct for autocorrelation with GLS, my $$\beta$$ estimates don't get any better. When my problem is slightly different and $$y_t= \sum_{i=1}^p \beta_i x_t^i + u_t,$$ such that $$y$$ is only a function of different $$x$$'s at the concurrent time, GLS does amazing, while least squares struggles.

Why does GLS perform so poorly in the first example and how do I get it to perform better (I only care about parameter estimation)?

I coded my two examples in R:

n_sim <- 15
LS_R2 <- rep(0, n_sim)
GLS_R2 <- rep(0, n_sim)
ar_noise <- 0.7
for(sim in 1:n_sim){
n=5000
set.seed(sim)
x_vec <- arima.sim(list(order=c(1,0,0), ar=c(0.5)), n=n, sd=1)
noise <- arima.sim(list(order=c(1,0,0), ar=c(ar_noise)), n=n, sd=1)
y <- rep(0,n)
p <- 150
true_beta <- seq(from=5, to=0, length.out=p) + rnorm(p,sd=0.1)
for(i in p:length(y)){
y[i] <- y[i] + sum(x_vec[i:(i-p+1)]*true_beta)
}
y <- y + noise*var(y)/var(noise)/50

# try simple least squares (should be very inefficient)
Xmat <- matrix(0, ncol = p, nrow = n)
for(i in 1:ncol(Xmat)){
Xmat[,i] <- c(rep(0,i-1), x_vec[1:(length(x_vec)-i+1)])
}

XXProd <- crossprod(Xmat)
XyProd <- crossprod(Xmat, matrix(y, ncol = 1))

est_beta <- as.numeric(solve(XXProd, XyProd, tol=0))
LS_R2[sim] <- 1- sum((true_beta-est_beta)^2)/sum((true_beta -
mean(true_beta))^2)

# try generalized least squares (should be optimal and efficient)
get_T <- function(ar, N){
i_vec <- c(1:N,2:N)
j_vec <- c(1:N,1:(N-1))
x_vec <- c(sqrt(1-ar^2), rep(1,N-1), rep(-ar,N-1))
sparseMatrix(i = i_vec, j = j_vec, x = x_vec, dims = c(N,N))
}

T1 <- get_T(ar_noise, length(y))

Xmat <- T1%*% Xmat
y <- as.numeric(T1%*% matrix(y, ncol=1))

XXProd <- crossprod(Xmat)
XyProd <- crossprod(Xmat, matrix(y, ncol = 1))

est_beta <- as.numeric(solve(XXProd, XyProd, tol=0))
GLS_R2[sim] <- 1- sum((true_beta-est_beta)^2)/sum((true_beta -
mean(true_beta))^2)
}
LS_R2
GLS_R2
mean(LS_R2)
mean(GLS_R2)


and

n_sim <- 15
LS_R2 <- rep(0,n_sim)
GLS_R2 <- rep(0,n_sim)
for(sim in 1:n_sim){
n=5000
set.seed(sim)
noise <- arima.sim(list(order=c(1,0,0), ar=c(ar_noise)), n=n, sd=1)
y <- rep(0,n)
p <- 150
Xmat <- matrix(0, ncol = p, nrow = n)
for(i in 1:ncol(Xmat)){
Xmat[,i] <- arima.sim(list(order=c(1,0,0), ar=c(0.5)), n=n, sd=1)
}
true_beta <- seq(from=5, to=0, length.out=p) + rnorm(p,sd=0.1)
for(i in p:length(y)){
y[i] <- y[i] + sum(Xmat[i,]*true_beta)
}
y <- y + noise*sd(y)/sd(noise)

# try simple least squares (should be very inefficient)
XXProd <- crossprod(Xmat)
XyProd <- crossprod(Xmat, matrix(y, ncol = 1))

est_beta <- as.numeric(solve(XXProd, XyProd, tol=0))
LS_R2[sim] <- 1- sum((true_beta-est_beta)^2)/sum((true_beta -
mean(true_beta))^2)

# try generalized least squares (should be optimal and efficient)
get_T <- function(ar,N){
i_vec <- c(1:N,2:N)
j_vec <- c(1:N,1:(N-1))
x_vec <- c(sqrt(1-ar^2), rep(1,N-1), rep(-ar,N-1))
sparseMatrix(i = i_vec, j = j_vec, x = x_vec, dims = c(N,N))
}

T1 <- get_T(ar_noise, length(y))

Xmat <- T1 %*% Xmat
y <- as.numeric(T1%*% matrix(y,ncol=1))

XXProd <- crossprod(Xmat)
XyProd <- crossprod(Xmat, matrix(y,ncol = 1))

est_beta <- as.numeric(solve(XXProd, XyProd,tol=0))
GLS_R2[sim] <- 1- sum((true_beta-est_beta)^2)/sum((true_beta -
mean(true_beta))^2)

}
LS_R2
GLS_R2
mean(LS_R2)
mean(GLS_R2)

• In the first one, you are using the same $x$ at different lags so the $x$ have a large chance of being correlated because they are the same variable $x_t$ at different points adjacent in time. In the second one, you are using different $x's$ so they are probably much less correlated. Note that you need to be careful between the correlation of the independent variables with themselves. That's not the same as the autocorrelation of the residuals ? So, in the first case, you may have a case of multi-collinearity in which case, one would not expect GLS to help. Commented Jul 7 at 17:19
• but in both cases I have autocorrelated residuals. Commented Jul 7 at 20:03
• It sounds interesting but I don't have much time, otherwise I'd try to look at it and help. But, for the issue with the first one, check out the koyck distributed lag if you haven't yet. That ( most decent econometrics texts will have a section covering distributed lags of which that is a specific case ) will give you some ideas for how to possibly fix that issue. That's as much as I can guess ( if it doesn't fix it, write back again ) but that discussion will lead to Almon lag etc. Commented Jul 8 at 1:34
• For the second one, I've never done that and I'm not sure that I follow it. You have $p$ different independent variables that are all measured at the same time $t$ and the response is also measured at this time. Then, the data set is whatever total number of observations of that ? If that's right, then I guess you'd need to make sure that you don't need lags of those independent variables. It sounds like the response currently just turns on ( do you need different shocks for each of the p independents ? ) and then immediately turns off which seems odd. Commented Jul 8 at 1:36
• the last sentence of immediately above should read " turns on at time period $t$ and turns off in the next time period $(t+1)$ which seems odd". Commented Jul 8 at 1:49