Software 4.7.1 (A Monte Carlo study of least squares bias). A Monte Carlo experiment answers questions that theory answers only approximately: generate many samples from a known model, estimate on each, and look at the distribution of the estimates. Here \(Y_t = 1 + 0.75Y_{t-1} + u_t\) with \(u_t \sim N(0, 1)\), 2 000 samples for each length \(T\):

> set.seed(2026)
> ar1_ols <- function(T, beta = 0.75, alpha = 1) {
+   y <- numeric(T + 50); y[1] <- alpha / (1 - beta)
+   for (t in 2:(T + 50)) y[t] <- alpha + beta * y[t - 1] + rnorm(1)
+   y <- y[-(1:50)]                     # discard a burn-in
+   coef(lm(y[-1] ~ y[-T]))[2] }
> for (T in c(20, 50, 200, 1000)) {
+   b <- replicate(2000, ar1_ols(T))
+   cat(sprintf("T = %4d  mean b = %.4f  approx bias = %.4f  sd = %.4f\n",
+       T, mean(b), -(1 + 3 * 0.75) / T, sd(b))) }
T =   20  mean b = 0.5796  approx bias = -0.1625  sd = 0.1949
T =   50  mean b = 0.6778  approx bias = -0.0650  sd = 0.1143
T =  200  mean b = 0.7326  approx bias = -0.0163  sd = 0.0500
T = 1000  mean b = 0.7479  approx bias = -0.0032  sd = 0.0209

Reading the output. With twenty observations the average estimate is \(0.580\), a bias of \(-0.170\), close to the approximation \(-0.1625\). The bias falls roughly in proportion to \(1/T\), and by \(T = 1000\) it is negligible: this is consistency at work. The standard deviation falls like \(1/\sqrt T\). The asymptotic value \(\sqrt {(1 - 0.75^2)/T}\) is \(0.047\) at \(T = 200\) and \(0.021\) at \(T = 1000\), close to the simulated \(0.050\) and \(0.021\). The first 50 values of each series are discarded so that the process has forgotten its starting value. That “burn-in” is standard practice whenever a stationary process is simulated.

Stata: simulate b = r(b), reps(2000): myprog with a program that generates one series and regresses it. EViews: a program loop with for and nrnd, storing @coefs(2) in a vector.

Software 4.7.2 (Least squares, conditional and exact maximum likelihood). One simulated series of 100 observations from an AR(1) with \(\beta = 0.75\) and mean \(4\), estimated three ways:

> set.seed(7)
> y <- arima.sim(list(ar = 0.75), n = 100) + 4
> coef(lm(y[-1] ~ y[-100]))
(Intercept)     y[-100]
  1.1509721   0.7310427
> fit_css <- arima(y, order = c(1, 0, 0), method = "CSS")
> fit_ml  <- arima(y, order = c(1, 0, 0), method = "ML")
> rbind(CSS = coef(fit_css), ML = coef(fit_ml))
          ar1 intercept
CSS 0.7310431  4.279385
ML  0.7548207  4.461537

Reading the output. Conditional maximum likelihood (CSS, “conditional sum of squares”) reproduces least squares to six decimals, as Theorem 4.6.1(i) says it must. Note that arima reports the process mean, labelled intercept, not \(\alpha \). The two are linked by \(\mu = \alpha /(1 - \beta )\), and indeed \(1.1510/(1 - 0.7310) = 4.279\). Exact maximum likelihood, which also uses the distribution of the first observation, gives the somewhat different \(\hat \beta = 0.755\). With \(T = 100\) the gap is well within one standard error, about \(0.07\).

Stata: arima y, ar(1) (exact ML by default; condition for conditional). EViews: ls y c ar(1), which estimates by non-linear least squares; the ARMA options choose conditional or exact ML.