1
$\begingroup$

Let’s start loading fpp3 package (https://github.com/robjhyndman/fpp3-package) and the US consumption expenditure dataset (https://rdrr.io/cran/fpp3/man/us_change.html)

library(fpp3)

us_change <- readr::read_csv("https://otexts.com/fpp3/extrafiles/us_change.csv") %>%
  mutate(Time = yearquarter(Time)) %>%
  as_tsibble(index = Time)

Suppose we want to forecast changes in Consumption based on changes in Income, so we use Income as predictor.

Let's start with a simple linear regression model

fit_lm <- us_change %>% model(TSLM(Consumption ~ Income))

Model report:

> report(fit_lm)
Series: Consumption 
Model: TSLM 

Residuals:
     Min       1Q   Median       3Q      Max 
-2.40845 -0.31816  0.02558  0.29978  1.45157 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.54510    0.05569   9.789  < 2e-16 ***
Income       0.28060    0.04744   5.915 1.58e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.6026 on 185 degrees of freedom
Multiple R-squared: 0.159,  Adjusted R-squared: 0.1545
F-statistic: 34.98 on 1 and 185 DF, p-value: 1.5774e-08

A sample of residuals:

> residuals(fit_lm)
# A tsibble: 187 x 3 [1Q]
# Key:       .model [1]
   .model                        Time  .resid
   <chr>                        <qtr>   <dbl>
 1 TSLM(Consumption ~ Income) 1970 Q1 -0.202 
 2 TSLM(Consumption ~ Income) 1970 Q2 -0.413 
 3 TSLM(Consumption ~ Income) 1970 Q3 -0.104 
 4 TSLM(Consumption ~ Income) 1970 Q4 -0.748 
 5 TSLM(Consumption ~ Income) 1971 Q1  0.795 
 6 TSLM(Consumption ~ Income) 1971 Q2 -0.0392
 7 TSLM(Consumption ~ Income) 1971 Q3  0.100 
 8 TSLM(Consumption ~ Income) 1971 Q4  0.778 
 9 TSLM(Consumption ~ Income) 1972 Q1  0.640 
10 TSLM(Consumption ~ Income) 1972 Q2  1.06  
# ... with 177 more rows

Now we suppose that the errors from the regression contain autocorrelation, so we use a regression model with ARIMA errors

fit_reg_arima_errs <- us_change %>% model(ARIMA(Consumption ~ Income))

Model report:

> report(fit_reg_arima_errs)
Series: Consumption 
Model: LM w/ ARIMA(1,0,2) errors 

Coefficients:
         ar1      ma1     ma2  Income  intercept
      0.6922  -0.5758  0.1984  0.2028     0.5990
s.e.  0.1159   0.1301  0.0756  0.0461     0.0884

sigma^2 estimated as 0.3219:  log likelihood=-156.95
AIC=325.91   AICc=326.37   BIC=345.29

In this case we have two types of residuals: the residuals from the regression model and the residuals from the ARIMA model.

Sample of regression residuals:

> regression_errors = residuals(fit_reg_arima_errs, type="regression") %>% print()
# A tsibble: 187 x 3 [1Q]
# Key:       .model [1]
   .model                         Time .resid
   <chr>                         <qtr>  <dbl>
 1 ARIMA(Consumption ~ Income) 1970 Q1  0.616
 2 ARIMA(Consumption ~ Income) 1970 Q2  0.460
 3 ARIMA(Consumption ~ Income) 1970 Q3  0.877
 4 ARIMA(Consumption ~ Income) 1970 Q4 -0.274
 5 ARIMA(Consumption ~ Income) 1971 Q1  1.90 
 6 ARIMA(Consumption ~ Income) 1971 Q2  0.912
 7 ARIMA(Consumption ~ Income) 1971 Q3  0.795
 8 ARIMA(Consumption ~ Income) 1971 Q4  1.65 
 9 ARIMA(Consumption ~ Income) 1972 Q1  1.31 
10 ARIMA(Consumption ~ Income) 1972 Q2  1.89 
# ... with 177 more rows

Sample of ARIMA residuals:

ARIMA_errors = residuals(fit_reg_arima_errs, type="innovation") %>% print()
# A tsibble: 187 x 3 [1Q]
# Key:       .model [1]
   .model                         Time  .resid
   <chr>                         <qtr>   <dbl>
 1 ARIMA(Consumption ~ Income) 1970 Q1 -0.167 
 2 ARIMA(Consumption ~ Income) 1970 Q2 -0.320 
 3 ARIMA(Consumption ~ Income) 1970 Q3  0.0720
 4 ARIMA(Consumption ~ Income) 1970 Q4 -0.694 
 5 ARIMA(Consumption ~ Income) 1971 Q1  1.05  
 6 ARIMA(Consumption ~ Income) 1971 Q2  0.142 
 7 ARIMA(Consumption ~ Income) 1971 Q3 -0.0525
 8 ARIMA(Consumption ~ Income) 1971 Q4  0.695 
 9 ARIMA(Consumption ~ Income) 1972 Q1  0.469 
10 ARIMA(Consumption ~ Income) 1972 Q2  0.788 
# ... with 177 more rows

Shouldn't the regression residuals from this latter model be identical (or at least similar) to the residuals from the first linear regression model?

And why do the regression residuals from the regression with ARIMA errors coincide with the values ​​of the response variable (Consumption)?

> us_change[,c("Time","Consumption")]
# A tsibble: 187 x 2 [1Q]
      Time Consumption
     <qtr>       <dbl>
 1 1970 Q1       0.616
 2 1970 Q2       0.460
 3 1970 Q3       0.877
 4 1970 Q4      -0.274
 5 1971 Q1       1.90 
 6 1971 Q2       0.912
 7 1971 Q3       0.795
 8 1971 Q4       1.65 
 9 1972 Q1       1.31 
10 1972 Q2       1.89 
# ... with 177 more rows

What am I missing?

The sample code is taken from "Forecasting: Principles and Practice. / Hyndman, Robin John; Athanasopoulos, George." (https://otexts.com/fpp3/)

Link for the dataset: https://otexts.com/fpp3/extrafiles/us_change.csv

$\endgroup$
  • 1
    $\begingroup$ One should model the actual observed Y and X and then predict a future Y using a future X AND then compute/impute the percent change objective that you are after. I don't have routine/easy access to the data that you are using so you can either email it to me and/or actually include it specifically in your post as a column oriented csv file and time available I will try and duplicate your results and deal with your questions in a subsequent response. $\endgroup$ – IrishStat Dec 19 '19 at 15:27
  • $\begingroup$ Thank you for your interest! I edited the post including a link (at the bottom of the post) to the .csv dataset used as sample data. A clarification, my goal now is only residuals analysis, I'm not actually forecasting $\endgroup$ – Raffaele Giannella Dec 19 '19 at 15:47
4
$\begingroup$
  1. A regression model with ARIMA errors is estimated using maximum likelihood estimation. It will give different results from estimating a regression model using least squares estimation because the optimization criterion is different.

  2. There was a bug in the fable package for computing regression residuals. Now fixed in the github version (https://github.com/tidyverts/fable).

$\endgroup$
  • $\begingroup$ Thank you Prof. Hyndman, I understood the difference between the optimization criterions and I confirm that using the github version of the Fable package my code works as I expected. $\endgroup$ – Raffaele Giannella Dec 20 '19 at 9:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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