how to get residuals from an autoregressive model by hand

Question

I would like to know how to calculate residuals from a fitted AR model by hand. Now I can use stats::ar to fit an autoregresive model and and get the residuals by from $resid of the fitted model, however I would like to be able to calculate these myself given the estimated AR parameters. So I would like to have a function that takes the estimated AR parameters from the ar model and a time-series and returns the residuals or a time-series with the auto-correlation removed.

My ultimate goal is to be able to estimate AR coefficients on one dataset and apply them to another, or to try to remove slightly different AR coefficients from the same time-series.

So for example here I would lieke to be able to reproduce the same ar_model$resid data by applying the estimated coefs to the original multivariate usconsumption timeseries.

> library(fpp)
> ar_model <- ar(usconsumption, aic=FALSE, order.max=1)
> ar_model$ar
, , consumption

  consumption    income
1   0.3088803 0.5627686

, , income

  consumption     income
1  0.08269438 -0.2310507

> head(ar_model$resid)
         consumption     income
1970 Q1          NA         NA
1970 Q2  -0.2363646  1.0249511
1970 Q3   0.1294457  1.0084068
1970 Q4  -1.1150108 -0.9913129
1971 Q1   1.5423841  1.5613124
1971 Q2  -0.2947244  0.3983440

I tried to use timsac::mfilter function, but I wasn't able to get anything close to the original ar residuals.

Answer 1

I wrote this code to try to answer your question. It will multiply the AR coefficients by the lagged values according to this formula from the ar help page.

Unfortunately, the result is not exactly correct :( My calculation of the predicted values is not the same as fitted(model).

Although this does not answer your question, I will post it here in case other users are able identify what's wrong wit this code, or use it as a template for a full solution.

formula:
x t ​ −μ=a 1 ​ (x t−1 ​ −μ)+⋯+a p ​ (x t−p ​ −μ)+e t ​

library(fpp, quietly = TRUE)

# fit
model <- ar(usconsumption, aic = FALSE, order.max = 1)

# initialize
pred <- matrix(nrow = nrow(usconsumption), ncol = 2)
colnames(pred) <- colnames(usconsumption)

# predict
for (i in 2:nrow(usconsumption)) {
  pred[i, "consumption"] <- model$x.mean["consumption"] + sum((usconsumption[i-1, ] - model$x.mean) * model$ar[, ,"consumption"])
  pred[i, "income"] <- model$x.mean["income"] + sum((usconsumption[i-1, ] - model$x.mean) * model$ar[, ,"income"])
}

# compare
head(fitted(model))
#>         getResponse(object).consumption getResponse(object).income
#> 1970 Q1                              NA                         NA
#> 1970 Q2                       0.6912944                  0.7115085
#> 1970 Q3                       0.7452274                  0.3364742
#> 1970 Q4                       0.8424964                  0.6631670
#> 1971 Q1                       0.3498029                  0.4041199
#> 1971 Q2                       1.2081025                  1.0924131

head(pred)
#>      consumption    income
#> [1,]          NA        NA
#> [2,]   0.5760684 0.7801837
#> [3,]   1.2252548 0.4806877
#> [4,]   1.1345370 0.6058726
#> [5,]  -0.1613336 0.8975607
#> [6,]   1.7980539 0.5466365

^{Created on 2023-05-22 with reprex v2.0.2}