2
$\begingroup$

I am conducting a Monte Carlo simulation to assess how OLS homoskedastic standard errors (i.e., based on $\sigma^2(X^\prime X)^{-1}$) change as serial correlation changes in a bivariate linear regression. Theoretically, OLS standard errors should be too large (relative to the true standard deviation of the estimate) when the disturbances are negatively correlated and too small when they are positively correlated.

However, this is not what my simulation is producing. What I am finding is shown in Figure 1, where OLS standard errors are generally too small with negative serial correlation and too large with positive serial correlation- the exact opposite of what is theoretically expected.

My simplified simulation is given below. Serial correlation can be varied by passing a different number to p in the set of {-0.99, 0.99}. The number of simulations and observations have been kept low for demonstration purposes. The results object contains the output of the simulation and some post-simulation calculations, like bias. In my full simulation, I also calculated Newey-West autocorrelation consistent standard errors and they were similar to the OLS standard errors, which troubles me further. In addition, I have considered other ways to induce serial correlation, such as using stats::filter and manually creating an autocorrelated series with another for loop.

What am I doing wrong here? Any help would be greatly appreciated. Again, the minimal goal of the simulation is to show how OLS standard errors are overconfident (i.e., smaller than the standard deviation of the estimate) as serial correlation increases and is positive. For background, this is for teaching purposes and not for a homework assignment.

set.seed(1)
M=10
N=100
p=0.99
data<-data.frame(
  x=rnorm(N)
)
results<-data.frame(
  b1_hat=rep(NA,M),
  se=rep(NA,M),
  sd=rep(NA,M),
  confidence=rep(NA,M),
  bias=rep(NA,M)
)
for (m in 1:M) {
  e<-arima.sim(list(ar=p),N)
  for (i in 1:N) {
    data$y[i]<-1+2*data$x[i]+e[i]
  }
  results[m,1]<-summary(lm(y~x,data))$coefficients[2]
  results[m,2]<-summary(lm(y~x,data))$coefficients[4]
}
results[,3]<-sqrt(sum((results$b1_hat-mean(results$b1_hat))^2)/length(results$b1_hat))
results[,4]<-results$sd/mean(results$se)
results[,5]<-mean(results$b1_hat)-2
$\endgroup$
3
  • $\begingroup$ (1) M is far too small to estimate standard errors accurately. (2) Your figure doesn't look like what you describe: it shows SEs increasing as one moves away from 0 in either direction. This will become more evident when you increase M to, say, 100 or more. (3) By using the built-in vector operations in R you will clarify your code and make it more efficient and reliable. (4) Your formulas in the last lines don't accomplish what you claim they do. Just plot the standard deviation of the estimates against p. $\endgroup$
    – whuber
    Commented Nov 13, 2023 at 20:04
  • $\begingroup$ Compare the output of this little study: M <- 100; data <- data.frame(x = seq(-3, 3, length.out = N)); data$y0 <- 1 + 2 * data$x; rho <- seq(-0.99, 0.99, length.out = 21); se <- sapply(rho, \(p) { se <- replicate(M, { data$y <- data$y0 + arima.sim(list(ar = p), N); coefficients(lm(y ~ x, data))["x"] }); sd(se) }); plot(rho, se, type = "l") $\endgroup$
    – whuber
    Commented Nov 13, 2023 at 20:09
  • $\begingroup$ Thank you for responding. You are correct that ideally more trials should be run. When I run the same simulation with 500 trials (M) and 1000 observations per model (N), I still find that the standard errors are larger than the standard deviations of the estimate when there is positive serial correlation. It is unclear what is incorrect about my calculations but I would love to hear more. Also, I do not understand why plotting standard deviations of the estimates against p is instructive here. Shouldn't the standard deviations of the estimates should be stable despite serial correlation? $\endgroup$ Commented Nov 13, 2023 at 23:16

1 Answer 1

1
$\begingroup$

I figured out the problem. As Baltagi (2008) states in eq. 5.34 (p. 111), the estimated variance of the point estimate when serial correlation exists is the original estimator (i.e., s^2(X'X)^{-1}) plus the product of the serial correlation, autocorrelation of the sample data, and the variance of the residuals. Since I generated my sample data as white noise, the second term in eq. 5.34 would be zero. Thus, adding serial correlation by itself was insufficient to produce the theoretically expected results.enter image description here

Once making the sample data autocorrelated as well (AR(0.5)), the results of my full simulation comports to theoretical expecations. OLS standard errors tend to overestimate the true standard deviation of the point estimate when negative serial correlation is present and underestimate it when positive serial correlation is present. The results can be seen in the plot below. Also, Newey-West autocorrelation consistent standard errors now tend to correct the issue but not well for extreme values of AR(p).

The full code I used for my simulation is given below.

mc<-function(p){
  library(sandwich)
  set.seed(1)
  M=500
  N=1000
  data<-data.frame(
    x=arima.sim(list(ar=0.5),N)
  )
  results<-data.frame(
    b1_hat=rep(NA,M),
    se_OLS=rep(NA,M),
    se_NW=rep(NA,M),
    sd=rep(NA,M),
    confidence=rep(NA,M),
    bias=rep(NA,M),
    rho=rep(p,M)
  )
  for (m in 1:M) {
    e<-arima.sim(list(ar=p),N)
    for (i in 1:N) {
      data$y[i]<-1+2*data$x[i]+e[i]
    }
    results[m,1]<-summary(lm(y~x,data))$coefficients[2]
    results[m,2]<-summary(lm(y~x,data))$coefficients[4]
    results[m,3]<-sqrt(diag(NeweyWest(lm(y~x,data),lag=(N^(1/4)))))[2]
  }
  results[,4]<-sqrt(sum((results$b1_hat-mean(results$b1_hat))^2)/length(results$b1_hat))
  results[,5]<-results$sd/mean(results$se_OLS)
  results[,6]<-mean(results$b1_hat)-2
  return(results)
}
a<-mc(-0.99)
b<-mc(-0.75)
c<-mc(-0.25)
d<-mc(-0.5)
e<-mc(0)
f<-mc(0.25)
g<-mc(0.5)
h<-mc(0.75)
i<-mc(0.99)
library(tidyverse)
mcdata<-full_join(a,b) %>%
  full_join(.,c) %>%
  full_join(.,d) %>%
  full_join(.,e) %>%
  full_join(.,f) %>%
  full_join(.,g) %>%
  full_join(.,h) %>%
  full_join(.,i)
plot<-ggplot(mcdata,aes(rho,se_OLS))+geom_point()
plot<-plot+geom_line(aes(rho,sd))
plot<-plot+labs(title="The effect of serial correlation on standard errors.",
                subtitle="500 simulated models each with 1000 observations.\nPoints represent the standard errors of the slope estimate from each model.\nThe line represents the standard deviation of the slope estimate for each value of AR(p).",
                x="Order of AR(p) disturbances",
                y="OLS standard errors")

Baltagi, Badi. 2008. Econometrics. Springer. DOI 10.1007/978-3-540-76516-5.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.