Naive, Basic, Dynamic, and Advanced Forecasting Models
Sith
2023-11-16
Libraries and Data
require(fpp3)## Loading required package: fpp3## Warning: package 'fpp3' was built under R version 4.3.2## ── Attaching packages ────────────────────────────────────────────── fpp3 0.5 ──## ✔ tibble 3.2.1 ✔ tsibble 1.1.3
## ✔ dplyr 1.1.4 ✔ tsibbledata 0.4.1
## ✔ tidyr 1.3.0 ✔ feasts 0.3.1
## ✔ lubridate 1.9.3 ✔ fable 0.3.3
## ✔ ggplot2 3.4.4 ✔ fabletools 0.3.4## Warning: package 'tibble' was built under R version 4.3.2## Warning: package 'dplyr' was built under R version 4.3.2## Warning: package 'tidyr' was built under R version 4.3.2## Warning: package 'lubridate' was built under R version 4.3.2## Warning: package 'ggplot2' was built under R version 4.3.2## Warning: package 'tsibble' was built under R version 4.3.2## Warning: package 'tsibbledata' was built under R version 4.3.2## Warning: package 'feasts' was built under R version 4.3.2## Warning: package 'fabletools' was built under R version 4.3.2## Warning: package 'fable' was built under R version 4.3.2## ── Conflicts ───────────────────────────────────────────────── fpp3_conflicts ──
## ✖ lubridate::date() masks base::date()
## ✖ dplyr::filter() masks stats::filter()
## ✖ tsibble::intersect() masks base::intersect()
## ✖ tsibble::interval() masks lubridate::interval()
## ✖ dplyr::lag() masks stats::lag()
## ✖ tsibble::setdiff() masks base::setdiff()
## ✖ tsibble::union() masks base::union()require(ggplot2)
require(kableExtra)## Loading required package: kableExtra## Warning: package 'kableExtra' was built under R version 4.3.2## Error: package or namespace load failed for 'kableExtra':
## .onLoad failed in loadNamespace() for 'kableExtra', details:
## call: !is.null(rmarkdown::metadata$output) && rmarkdown::metadata$output %in%
## error: 'length = 2' in coercion to 'logical(1)'require(latex2exp)## Loading required package: latex2exp## Warning: package 'latex2exp' was built under R version 4.3.2require(readxl)## Loading required package: readxl## Warning: package 'readxl' was built under R version 4.3.2require(tidyverse)## Loading required package: tidyverse## Warning: package 'tidyverse' was built under R version 4.3.2## Warning: package 'readr' was built under R version 4.3.2## Warning: package 'purrr' was built under R version 4.3.2## Warning: package 'stringr' was built under R version 4.3.2## Warning: package 'forcats' was built under R version 4.3.2## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats 1.0.0 ✔ readr 2.1.4
## ✔ purrr 1.0.2 ✔ stringr 1.5.1## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ tsibble::interval() masks lubridate::interval()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorsrequire(magrittr)## Loading required package: magrittr## Warning: package 'magrittr' was built under R version 4.3.2##
## Attaching package: 'magrittr'
##
## The following object is masked from 'package:purrr':
##
## set_names
##
## The following object is masked from 'package:tidyr':
##
## extractrequire(fitdistrplus)## Loading required package: fitdistrplus## Warning: package 'fitdistrplus' was built under R version 4.3.3## Loading required package: MASS
##
## Attaching package: 'MASS'
##
## The following object is masked from 'package:dplyr':
##
## select
##
## Loading required package: survival## Warning: package 'survival' was built under R version 4.3.2library(fredr)## Warning: package 'fredr' was built under R version 4.3.2require(MASS)
require(psych)## Loading required package: psych## Warning: package 'psych' was built under R version 4.3.2##
## Attaching package: 'psych'
##
## The following objects are masked from 'package:ggplot2':
##
## %+%, alpharequire(ResourceSelection)## Loading required package: ResourceSelection## Warning: package 'ResourceSelection' was built under R version 4.3.2## ResourceSelection 0.3-6 2023-06-27require(eia)## Loading required package: eia## Warning: package 'eia' was built under R version 4.3.2## EIA_KEY not found. See `vignette("api", "eia")` for key storage options.Keys
## Key stored successfully in package environment.Read & Edit Data
total=fredr(
series_id = "USACPALTT01CTGYM",
observation_start = as.Date("2010-01-01"),
observation_end = as.Date("2023-12-31")
)
treas=fredr(
series_id = "TB3MS",
observation_start = as.Date("2010-01-01"),
observation_end = as.Date("2023-12-31")
)
total$TBill=treas$value
total$Inflation=total$value
total$value=NULLBuild Tsibble
total=total%>%mutate(YearMonth=yearmonth(date))%>%as_tsibble(index=YearMonth)
total=total[,5:7]
temp=reshape2::melt(total, id.var='YearMonth')
ggplot(data=temp, aes(x=YearMonth, y=value, col=variable))+geom_line()+labs(col='Variable')
total$ElectionSpline=c(rep(0,12), rep(1,12), rep(0,48), rep(1,12), rep(0,48), rep(1,12), rep(0,24))
# geom_line(aes(x=YearMonth, y=Treas))
Regression with ARIMA Errors
Time Series and Train / Test Set
train=total[1:156,]
test=total[157:nrow(total),]
# Build an MMS variable
mymin=min(train$Inflation)
mymax=max(train$Inflation)
myrange=mymax-mymin
train$newInf=(train$Inflation-mymin)/myrange
test$newInf=(test$Inflation-mymin)/myrange
# Build an offset for TBill
temp=c(total$TBill[1:12], total$TBill[1:156])
train$LagTBill=temp[1:156]
test$LagTBill=temp[157:nrow(total)]Time Series Plots
total%>%gg_season(Inflation)
total%>%gg_lag(Inflation,geom='point')
total%>%ACF(Inflation)%>%autoplot()
total%>%gg_season(TBill)
total%>%gg_lag(TBill,geom='point')
total%>%ACF(TBill)%>%autoplot()
Transformation
lambda=train%>%features(Inflation, features = guerrero) |> pull(lambda_guerrero)
train%>%autoplot(box_cox(Inflation, lambda[1])) +
labs(y = "",title = latex2exp::TeX(paste0(
"Transformed Inflation with $\\lambda$ = ", round(lambda[1],2))))
lambda=train%>%features(TBill, features = guerrero) |> pull(lambda_guerrero)
train%>%autoplot(box_cox(TBill, lambda[1])) +
labs(y = "",title = latex2exp::TeX(paste0(
"Transformed TBill with $\\lambda$ = ", round(lambda[1],2))))
# Decomposition
comp=total%>%model(stl=STL(Inflation))%>%components()
comp%>%as_tsibble() %>%autoplot(Inflation)+ geom_line(aes(y=trend), colour = "red")+
geom_line(aes(y=season_adjust), colour = "blue")
comp%>%autoplot()
comp=total%>%model(stl=STL(TBill))%>%components()
comp%>%as_tsibble() %>%autoplot(TBill)+ geom_line(aes(y=trend), colour = "red")+
geom_line(aes(y=season_adjust), colour = "blue")
comp%>%autoplot()
Models
#Naive Models
m1=train |>model(Mean_Model=MEAN(Inflation))
m2=train |>model(Drift_Model=RW(Inflation~drift()))
m3=train |>model(Naive_Model=NAIVE(Inflation))
m4=train |>model(Seasonal_Naive_Model=SNAIVE(Inflation))
#Basic Models
m5=train%>%model(TS_Model = TSLM(Inflation))
m6=train |>model(ETS_Model=ETS(Inflation))
m7=train |>model(ARIMA_Model=ARIMA(Inflation))
#Dynamic Models
m8=train |>model(Fourier_Model=TSLM(Inflation~ElectionSpline+trend(knots=yearmonth('2022 Jun'))+fourier(K=2)))
m9=train |>model(Exponential_Model=TSLM(log(Inflation+abs(mymin)+.001)~ElectionSpline+trend(knots=yearmonth('2022 Jun'))+season()))
m10=train |>model(ARIMA(Inflation ~ ElectionSpline+LagTBill)) #ARIMA with Errors
#Dynamic Ensemble Models
m11=train|>model(Ensemble_Model=(ETS(Inflation)+ARIMA(Inflation~ElectionSpline+trend(knots = yearmonth('2022 Jun'))))/2)
m12=train|>model(Ensemble_Model_2=(ETS(Inflation)+TSLM(Inflation~trend(knots=yearmonth('2022 Jun'))+fourier(K=2))+ARIMA(Inflation))/3)
#Dynamic Neural Network Models
m13=train|>model(NeuralNet_Model=NNETAR(newInf~ElectionSpline+trend(knots=yearmonth('2022 Jun'))+season()))Example ARIMA + Regressor
fit <- train |>
model(ARIMA(Inflation ~ LagTBill))
report(fit)## Series: Inflation
## Model: LM w/ ARIMA(0,1,1)(0,0,1)[12] errors
##
## Coefficients:
## ma1 sma1 LagTBill intercept
## 0.2519 -0.2977 -0.4728 0.0416
## s.e. 0.0738 0.0942 0.2577 0.0265
##
## sigma^2 estimated as 0.1355: log likelihood=-63.58
## AIC=137.16 AICc=137.56 BIC=152.38 bind_rows(
`Regression residuals` =
as_tibble(residuals(fit, type = "regression")),
`ARIMA residuals` =
as_tibble(residuals(fit, type = "innovation")),
.id = "type"
) |>
mutate(
type = factor(type, levels=c(
"Regression residuals", "ARIMA residuals"))
) |>
ggplot(aes(x = YearMonth, y = .resid)) +
geom_line() +
facet_grid(vars(type))
fit |> gg_tsresiduals()
augment(fit) |>
features(.innov, ljung_box, dof = 4, lag = 12)## # A tibble: 1 × 3
## .model lb_stat lb_pvalue
## <chr> <dbl> <dbl>
## 1 ARIMA(Inflation ~ LagTBill) 5.68 0.683 forecast(fit, new_data = test) |>
autoplot(train) +
labs(y = "Percentage change")
Reports
for (i in 1:13){report(eval(parse(text=paste0('m',i))))}## Series: Inflation
## Model: MEAN
##
## Mean: 2.081
## sigma^2: 4.625
## Series: Inflation
## Model: RW w/ drift
##
## Drift: 0.0431 (se: 0.0324)
## sigma^2: 0.1626
## Series: Inflation
## Model: NAIVE
##
## sigma^2: 0.1626
## Series: Inflation
## Model: SNAIVE
##
## sigma^2: 3.2442
## Series: Inflation
## Model: TSLM
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.3078 -1.3407 -0.4304 0.1691 6.8874
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.0810 0.1722 12.09 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.151 on 155 degrees of freedom
## Series: Inflation
## Model: ETS(A,N,N)
## Smoothing parameters:
## alpha = 0.9998999
##
## Initial states:
## l[0]
## -0.1605545
##
## sigma^2: 0.1645
##
## AIC AICc BIC
## 510.2112 510.3691 519.3608
## Series: Inflation
## Model: ARIMA(0,1,2)(0,0,1)[12]
##
## Coefficients:
## ma1 ma2 sma1
## 0.3289 0.1268 -0.2655
## s.e. 0.0812 0.0795 0.0916
##
## sigma^2 estimated as 0.1373: log likelihood=-65.02
## AIC=138.05 AICc=138.31 BIC=150.22
## Series: Inflation
## Model: TSLM
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.27630 -1.03983 -0.01826 0.52417 4.91624
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) -0.290236 0.256757 -1.130
## ElectionSpline 0.209959 0.290824 0.722
## trend(knots = yearmonth("2022 Jun"))trend 0.028229 0.002841 9.937
## trend(knots = yearmonth("2022 Jun"))trend_150 0.793531 0.174087 4.558
## fourier(K = 2)C1_12 -0.029505 0.172567 -0.171
## fourier(K = 2)S1_12 0.137280 0.175245 0.783
## fourier(K = 2)C2_12 0.058285 0.172453 0.338
## fourier(K = 2)S2_12 0.017627 0.172694 0.102
## Pr(>|t|)
## (Intercept) 0.260
## ElectionSpline 0.471
## trend(knots = yearmonth("2022 Jun"))trend < 2e-16 ***
## trend(knots = yearmonth("2022 Jun"))trend_150 1.07e-05 ***
## fourier(K = 2)C1_12 0.864
## fourier(K = 2)S1_12 0.435
## fourier(K = 2)C2_12 0.736
## fourier(K = 2)S2_12 0.919
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.522 on 148 degrees of freedom
## Multiple R-squared: 0.5218, Adjusted R-squared: 0.4991
## F-statistic: 23.07 on 7 and 148 DF, p-value: < 2.22e-16
## Series: Inflation
## Model: TSLM
## Transformation: log(Inflation + `abs(mymin)` + 0.001)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.5978 -0.3955 0.1441 0.5826 1.7107
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) -0.909909 0.310344 -2.932
## ElectionSpline 0.152008 0.190350 0.799
## trend(knots = yearmonth("2022 Jun"))trend 0.014503 0.001859 7.800
## trend(knots = yearmonth("2022 Jun"))trend_150 0.118110 0.114353 1.033
## season()year2 0.080515 0.390722 0.206
## season()year3 0.092995 0.390735 0.238
## season()year4 -0.328262 0.390757 -0.840
## season()year5 0.027778 0.390788 0.071
## season()year6 0.148252 0.390828 0.379
## season()year7 0.157484 0.390901 0.403
## season()year8 0.166185 0.391157 0.425
## season()year9 0.087026 0.391594 0.222
## season()year10 0.143384 0.392212 0.366
## season()year11 0.174676 0.393010 0.444
## season()year12 0.285477 0.393987 0.725
## Pr(>|t|)
## (Intercept) 0.00393 **
## ElectionSpline 0.42588
## trend(knots = yearmonth("2022 Jun"))trend 1.25e-12 ***
## trend(knots = yearmonth("2022 Jun"))trend_150 0.30344
## season()year2 0.83703
## season()year3 0.81223
## season()year4 0.40229
## season()year5 0.94343
## season()year6 0.70501
## season()year7 0.68765
## season()year8 0.67159
## season()year9 0.82445
## season()year10 0.71523
## season()year11 0.65739
## season()year12 0.46991
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9961 on 141 degrees of freedom
## Multiple R-squared: 0.36, Adjusted R-squared: 0.2965
## F-statistic: 5.666 on 14 and 141 DF, p-value: 1.2515e-08
## Series: Inflation
## Model: LM w/ ARIMA(1,1,0)(0,0,1)[12] errors
##
## Coefficients:
## ar1 sma1 ElectionSpline LagTBill
## 0.3189 -0.2758 0.1506 -0.3845
## s.e. 0.0813 0.0903 0.1658 0.2624
##
## sigma^2 estimated as 0.1356: log likelihood=-63.62
## AIC=137.23 AICc=137.63 BIC=152.45
## Series: Inflation
## Model: COMBINATION
## Combination: Inflation * 0.5
##
## ============================
##
## Series: Inflation
## Model: COMBINATION
## Combination: Inflation + Inflation
##
## ==================================
##
## Series: Inflation
## Model: ETS(A,N,N)
## Smoothing parameters:
## alpha = 0.9998999
##
## Initial states:
## l[0]
## -0.1605545
##
## sigma^2: 0.1645
##
## AIC AICc BIC
## 510.2112 510.3691 519.3608
##
## Series: Inflation
## Model: LM w/ ARIMA(1,0,2)(0,0,1)[12] errors
##
## Coefficients:
## ar1 ma1 ma2 sma1 ElectionSpline
## 0.9728 0.3524 0.1355 -0.2501 0.1770
## s.e. 0.0228 0.0878 0.0822 0.0924 0.1625
## trend(knots = yearmonth("2022 Jun"))trend
## 0.0389
## s.e. 0.0131
## trend(knots = yearmonth("2022 Jun"))trend_150
## -0.3976
## s.e. 0.2211
##
## sigma^2 estimated as 0.1319: log likelihood=-61.88
## AIC=139.75 AICc=140.73 BIC=164.15
##
##
## Series: Inflation
## Model: COMBINATION
## Combination: Inflation * 0.333333333333333
##
## ==========================================
##
## Series: Inflation
## Model: COMBINATION
## Combination: Inflation + Inflation
##
## ==================================
##
## Series: Inflation
## Model: COMBINATION
## Combination: Inflation + Inflation
##
## ==================================
##
## Series: Inflation
## Model: ETS(A,N,N)
## Smoothing parameters:
## alpha = 0.9998999
##
## Initial states:
## l[0]
## -0.1605545
##
## sigma^2: 0.1645
##
## AIC AICc BIC
## 510.2112 510.3691 519.3608
##
## Series: Inflation
## Model: TSLM
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.32808 -1.05149 0.01303 0.47548 4.86270
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) -0.245028 0.248603 -0.986
## trend(knots = yearmonth("2022 Jun"))trend 0.028293 0.002835 9.980
## trend(knots = yearmonth("2022 Jun"))trend_150 0.780365 0.172851 4.515
## fourier(K = 2)C1_12 -0.028856 0.172287 -0.167
## fourier(K = 2)S1_12 0.134999 0.174935 0.772
## fourier(K = 2)C2_12 0.057842 0.172174 0.336
## fourier(K = 2)S2_12 0.016861 0.172413 0.098
## Pr(>|t|)
## (Intercept) 0.326
## trend(knots = yearmonth("2022 Jun"))trend < 2e-16 ***
## trend(knots = yearmonth("2022 Jun"))trend_150 1.28e-05 ***
## fourier(K = 2)C1_12 0.867
## fourier(K = 2)S1_12 0.442
## fourier(K = 2)C2_12 0.737
## fourier(K = 2)S2_12 0.922
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.52 on 149 degrees of freedom
## Multiple R-squared: 0.5201, Adjusted R-squared: 0.5007
## F-statistic: 26.91 on 6 and 149 DF, p-value: < 2.22e-16
##
##
## Series: Inflation
## Model: ARIMA(0,1,2)(0,0,1)[12]
##
## Coefficients:
## ma1 ma2 sma1
## 0.3289 0.1268 -0.2655
## s.e. 0.0812 0.0795 0.0916
##
## sigma^2 estimated as 0.1373: log likelihood=-65.02
## AIC=138.05 AICc=138.31 BIC=150.22
##
##
## Series: newInf
## Model: NNAR(2,1,9)[12]
##
## Average of 20 networks, each of which is
## a 17-9-1 network with 172 weights
## options were - linear output units
##
## sigma^2 estimated as 3.315e-05Plots
total$newInf=c(train$newInf,test$newInf)
myplot=function(model){
model|>forecast(test)|>autoplot(total)+geom_line(aes(y=.fitted,col=.model),
data=augment(model))+ggtitle(names(model))}
for (i in 1:13){print(myplot(eval(parse(text=paste0('m',i)))))} #NNETAR takes long to plot
## Warning: Removed 1 row containing missing values (`geom_line()`).
## Warning: Removed 1 row containing missing values (`geom_line()`).
## Warning: Removed 12 rows containing missing values (`geom_line()`).
## Warning: Removed 12 rows containing missing values (`geom_line()`).
Accuracy Metrics
tmp3=c(rep(0,12))
mymat=matrix(rep(0, 12*12), ncol=12)
for (i in 1:13){
tmp=eval(parse(text=paste0('m',i)))%>%forecast(test)
tmp2=tmp%>%accuracy(test)
tmp3=rbind(tmp2, tmp3)
}## Warning in rbind(deparse.level, ...): number of columns of result, 10, is not a
## multiple of vector length 12 of arg 2tmp3=tmp3[tmp3$RMSE>0,]
tmp3$MASE=tmp3$RMSSE=tmp3$.type=NULL
tmp3%>%arrange(MAPE)## # A tibble: 13 × 7
## .model ME RMSE MAE MPE MAPE ACF1
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Ensemble_Model -0.734 0.915 0.754 -21.3 21.6 0.530
## 2 ARIMA(Inflation ~ ElectionSpline … -0.948 1.11 0.978 -27.7 28.2 0.563
## 3 Mean_Model 2.01 2.28 2.01 45.8 45.8 0.674
## 4 TS_Model 2.01 2.28 2.01 45.8 45.8 0.674
## 5 ARIMA_Model -1.74 1.99 1.76 -51.0 51.3 0.690
## 6 NeuralNet_Model -0.226 0.244 0.226 -55.2 55.2 0.625
## 7 Naive_Model -2.43 2.67 2.43 -69.7 69.7 0.674
## 8 ETS_Model -2.43 2.67 2.43 -69.7 69.7 0.674
## 9 Drift_Model -2.71 2.97 2.71 -77.7 77.7 0.688
## 10 Seasonal_Naive_Model -3.89 4.10 3.89 -108. 108. 0.579
## 11 Ensemble_Model_2 -4.73 5.10 4.73 -134. 134. 0.721
## 12 Fourier_Model -10.1 10.8 10.1 -284. 284. 0.734
## 13 Exponential_Model -62.0 79.0 62.0 -1819. 1819. 0.659