ARIMA Homework
Jack Wright
9.1, 9.2, 9.3, 9.5, 9.6, 9.7, 9.8
9.1
Figure 9.32 shows the ACFs for 36 random numbers, 360 random numbers and 1,000 random numbers.
a. Explain the differences among these figures. Do they all indicate that the data are white noise?
## -- Attaching packages -------------------------------------------- fpp3 0.4.0 --
## v tibble 3.1.4 v tsibble 1.1.1
## v dplyr 1.0.7 v tsibbledata 0.4.0
## v tidyr 1.1.3 v feasts 0.2.2
## v lubridate 1.7.10 v fable 0.3.1
## v ggplot2 3.3.5
## -- Conflicts ------------------------------------------------- fpp3_conflicts --
## x lubridate::date() masks base::date()
## x dplyr::filter() masks stats::filter()
## x tsibble::intersect() masks base::intersect()
## x tsibble::interval() masks lubridate::interval()
## x dplyr::lag() masks stats::lag()
## x tsibble::setdiff() masks base::setdiff()
## x tsibble::union() masks base::union()
## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --
## v readr 2.0.1 v stringr 1.4.0
## v purrr 0.3.4 v forcats 0.5.1
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x lubridate::as.difftime() masks base::as.difftime()
## x lubridate::date() masks base::date()
## x dplyr::filter() masks stats::filter()
## x tsibble::intersect() masks lubridate::intersect(), base::intersect()
## x tsibble::interval() masks lubridate::interval()
## x dplyr::lag() masks stats::lag()
## x tsibble::setdiff() masks lubridate::setdiff(), base::setdiff()
## x tsibble::union() masks lubridate::union(), base::union()
These data are all white noise by definition (since we created the random data), the difference is the confidence intervals. The confidence intervals tell us at which point we deem the autocorrelation to be significant. This calculation is simply \(\pm 1.96/\sqrt{n}\), which we can confirm easily.
\[n=1500\\95\% ~confidence ~long = \pm1.96/\sqrt{1500}=.05\]
b. Why are the critical values at different distances from the mean of zero? Why are the autocorrelations different in each figure when they each refer to white noise?
As shown in part a., the critical values depend on the sample size. The autocorrelation values are different in each figure because we track the correlation between time t and time t-L. The longer the data is, the more that data converges to the true mean of the correlation. As you can see in the largest sample (plot 3) it is the furthest away from the critical values, because we are the most sure that this is white noise.
9.2
A classic example of a non-stationary series are stock prices. Plot the daily closing prices for Amazon stock (contained in gafa_stock), along with the ACF and PACF. Explain how each plot shows that the series is non-stationary and should be differenced.
gafa_stock%>%
filter(Symbol =='AMZN')%>%
autoplot(Close)+
labs(title= 'AMZN Closing Prices')
As we can see in this plot there is clear trend in the data. stationary data is trendless white noise.
Now lets look at the ACF for futher clues
amzn<-gafa_stock%>%
filter(Symbol =='AMZN')
amzn%>%
gg_tsdisplay(Close, plot_type = 'partial' )
## Warning: Provided data has an irregular interval, results should be treated with
## caution. Computing ACF by observation.
## Warning: Provided data has an irregular interval, results should be treated with
## caution. Computing ACF by observation.
As we can see, there is a significant trend to the ACF data. It looks exponential but the trend is so strong we don’t even really see the tail. This is a clear sign that differencing is required.
9.3
For the following series, find an appropriate Box-Cox transformation and order of differencing in order to obtain stationary data.
a. Turkish GDP from global_economy.
turkish<-global_economy%>%
filter(Country =='Turkey')%>%
select(Year,GDP)
turkish%>%
gg_tsdisplay(plot_type = 'partial')+
labs(title ='Untransfromed Turkish GDP')
## Plot variable not specified, automatically selected `y = GDP`
Both the plot of the raw data and the ACF plot show a clear trend in the data. Lets perform a first difference to the data and reixamine the ACF and PACF.
turkish%>%
gg_tsdisplay(difference(GDP), plot_type = 'partial')+
labs(title = 'Turkish GDP: First Difference, Untransformed')
## Warning: Removed 1 row(s) containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_point).
as we can see from the difference plot, the untransformed data’s variance is not constant throughout the series. Let’s see if a transformation helps.
library(forecast)
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
##
## Attaching package: 'forecast'
## The following object is masked from 'package:ggpubr':
##
## gghistogram
lambda<-forecast::BoxCox.lambda(turkish$GDP)
The suggested lambda is .15 . You might choose a log transformation for decipherability, but I prefer to use the actual lambda in this case.
transformed_turkish<-turkish%>%
mutate(GDP = forecast::BoxCox(GDP, lambda))
transformed_turkish%>%gg_tsdisplay(plot_type = 'partial')+
labs(title ='Transformed Turkish GDP')
## Plot variable not specified, automatically selected `y = GDP`
Now lets look at the first difference
transformed_turkish%>%
gg_tsdisplay(difference(GDP), plot_type = 'partial')+
labs(title = 'Turkish GDP: First Difference, Transformed, lambda =.15')
## Warning: Removed 1 row(s) containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_point).
As we can see, the variance is far more constant in the residuals, and both the ACF and PACF show no autocorrelation in the data.
Lets now use a unit root test on the transformed data to determine the number of differences required
transformed_turkish%>%
features(GDP, unitroot_ndiffs)
## # A tibble: 1 x 1
## ndiffs
## <int>
## 1 1
b.Accommodation takings in the state of Tasmania from aus_accommodation.
tasmania<-aus_accommodation%>%
filter(State =='Tasmania')%>%
select(Date, Takings)
tasmania%>%
gg_tsdisplay(plot_type = 'partial')+
labs(title ='Untransfromed Tasmanian Takings')
## Plot variable not specified, automatically selected `y = Takings`
Both the plot of the raw data and the ACF plot show a clear trend in the data as well as seasonality. Additionally the variance seems to increase with time as well. Since this question is only asking for a Box-Cox transformation, we won’t deal with the seasonality just yet.
tasmania%>%
gg_tsdisplay(difference(Takings), plot_type = 'partial')+
labs(title = 'Tasmanian Takings: First Difference, Untransformed')
## Warning: Removed 1 row(s) containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_point).
as we can see from the difference plot, the untransformed data’s variance is not constant throughout the series. Let’s see if a transformation helps.
library(forecast)
lambda<-BoxCox.lambda(tasmania$Takings)
lambda
## [1] -0.005076712
The suggested lambda is close to zero . I will choose a log transformation.
transformed_tasmania<-tasmania%>%
mutate(GDP = forecast::BoxCox(Takings, 1))
transformed_tasmania%>%gg_tsdisplay(plot_type = 'partial')+
labs(title ='Transformed Tasmanian Takings')
## Plot variable not specified, automatically selected `y = Takings`
Now lets look at the first difference
transformed_turkish%>%
gg_tsdisplay(difference(GDP), plot_type = 'partial')+
labs(title = 'Turkish GDP: First Difference, Transformed, lambda =.15')
## Warning: Removed 1 row(s) containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_point).
As we can see, the variance is far more constant in the residuals, and both the ACF and PACF show no autocorrelation in the data.
c. Monthly sales from souvenirs.
Since we have proven our process works, lets simplify this last example
lambda<-BoxCox.lambda(souvenirs$Sales)
souvenirs%>%
mutate(Sales = BoxCox(Sales, lambda))%>%
features(Sales, unitroot_ndiffs)
## # A tibble: 1 x 1
## ndiffs
## <int>
## 1 1
9.5
For your retail data (from Exercise 8 in Section 2.10), find the appropriate order of differencing (after transformation if necessary) to obtain stationary data.
Lets use our simplified process for this data as well.
set.seed(321)
myseries <- aus_retail %>%
filter(`Series ID` == sample(aus_retail$`Series ID`,1))
lambda<-BoxCox.lambda(myseries$Turnover)
print(lambda)
## [1] -0.2161443
myseries%>%
mutate(Turnover = BoxCox(Turnover, lambda))%>%
features(Turnover, unitroot_ndiffs)
## # A tibble: 1 x 3
## State Industry ndiffs
## <chr> <chr> <int>
## 1 Western Australia Newspaper and book retailing 1
9.6
Simulate and plot some data from simple ARIMA models.
*a. Use the following R code to generate data from an AR(1) model with
ϕ1 = 0.6 and \(σ^2\) = 1.
The process starts with
\(y_1\) = 0 *
y <- numeric(100)
e <- rnorm(100)
for(i in 2:100)
y[i] <- 0.6*y[i-1] + e[i]
sim <- tsibble(idx = seq_len(100), y = y, index = idx)
Produce a time plot for the series. How does the plot change as you change
\(ϕ_1\)?
Lets look at the time plot for the series first.
sim%>%autoplot()
## Plot variable not specified, automatically selected `.vars = y`
The data looks fairly random, without trend.
analyzing the equation presented, we see that y[i]’s value is dependent on a scalar \(\phi\) of the previously generated data plus a random error (ranging from about -2.5 to 2.5) In this context, it only makes sense to vary \(\phi\) from -1 to 1, because this clearly comes from an AR() model, in which the prediction for the next value is some weighted percentage of the previous values.
From our knowledge of AR() models:
1.When $\phi$ = 0, we will get white noise
2.When $\phi$ =1, we will get a random walk
3.When $\phi$ <0, $y_t$ will oscillate around the mean (in our case $\mu = 0$)
Write your own code to generate data from an MA(1) model with \(θ_1\)=0.6 and \(σ^2\)=1
The equation for a moving average model is as follows:
$$
y_t = c+_t + 1{t-1}+ 2{t-2}+++ q{t-q}
$$ Lets modify our function accordingly
y <- numeric(100)
sigma<-sqrt(1)
e <- rnorm(100, sigma)
for(i in 2:100)
y[i] <- 0.6*e[i-1] + e[i]
sim <- tsibble(idx = seq_len(100), y = y, index = idx)
sim%>%autoplot()
## Plot variable not specified, automatically selected `.vars = y`
Produce a time plot for the series. How does the plot change as you change
\(θ_1\)?
Since I am not aware of the rules of thumb for the actions of theta on a MA() model, I will do some simulation.
sim_MA<-function(theta){
y <- numeric(100)
sigma<-sqrt(1)
e <- rnorm(100, sigma)
e[1]<-0
for(i in 2:100)
y[i] <- theta*e[i-1] + e[i]
sim <- tsibble(idx = seq_len(100), y = y, index = idx)
}
get_mean<-function(sim,thet){
df<-NULL
for (i in 1:100){
temp<-sim(thet)
if(is.null(df)){
df<-as.vector(temp$y)
}else{
df<-cbind(df,as.vector(temp$y))
}
}
return (rowMeans(as.data.frame(df)))
}
create_phi_frame<-function(sim, theta){
data.frame(phi = as.character(theta),x =seq_len(100), test = get_mean(sim, theta))
}
sim_full_phi<-function(sim,min=-1,max=1){
phi_range<-seq(min,max,by=.1)
output_frame<-NULL
for (i in phi_range){
if(is.null(output_frame)){
output_frame = create_phi_frame(sim,theta = i)
}else{
output_frame<-rbind(output_frame,create_phi_frame(sim,theta=i))
}
}
return(output_frame)
}
dat<-sim_full_phi(sim_MA)
dat%>%group_by(phi)%>%summarize(mu = mean(test))%>%ggplot(aes(x=as.numeric(phi),y=mu))+geom_line()+ggtitle('Mean of MA() sim Data by theta value')
I simulated Data for each value of phi between -1 and 1 100 times and then took the mean of that data. It seems that as \(\theta\) increases, so does the average value of the data generated.
This makes analytical sense, because the \(\theta\) coefficient is acting as a x axis shift to the mean of the normally distributed error.
Generate data from an ARMA(1,1) model with \(ϕ_1\)=0.6 , \(θ_1\)=0.6 and \(σ^2\) = 1 .
I will use a similar technique of bootstrapping to get an average value for an ARMA(1,1)
sim_ARIMA11<-function(phi,theta){
y <- numeric(100)
sigma<-sqrt(1)
e <- rnorm(100, sigma)
e[1]<-0
for(i in 2:100)
y[i] <- phi*y[i-1]+theta*e[i-1]+e[i]
sim <- tsibble(idx = seq_len(100), y = y, index = idx)
}
get_mean2<-function(sim,ph,thet){
df<-NULL
for (i in 1:100){
temp<-sim(ph,thet)
if(is.null(df)){
df<-as.vector(temp$y)
}else{
df<-cbind(df,as.vector(temp$y))
}
}
return (rowMeans(as.data.frame(df)))
}
create_ARIMA_frame<-function(sim,phi, theta){
data.frame(model = 'ARIMA11',x =seq_len(100), test = get_mean2(sim, phi, theta))
}
Generate data from an AR(2) model with
ϕ 1 = − 0.8 ,
ϕ 2 = 0.3 and
σ 2 = 1 . (Note that these parameters will give a non-stationary series.)
I will repeat this strategy for the AR(2)
sim_AR2<-function(phi1,phi2){
y <- numeric(100)
sigma<-sqrt(1)
e <- rnorm(100, sigma)
e[1]<-0
for(i in 2:100)
y[i] <- phi1*y[i-1]+phi2*y[i-2]+e[i]
sim <- tsibble(idx = seq_len(100), y = y, index = idx)
}
get_mean3<-function(sim,ph1,ph2){
df<-NULL
for (i in 1:100){
temp<-sim(ph1,ph2)
if(is.null(df)){
df<-as.vector(temp$y)
}else{
df<-cbind(df,as.vector(temp$y))
}
}
return (rowMeans(as.data.frame(df)))
}
create_AR2_frame<-function(sim,phi1, phi2){
data.frame(model = 'AR2',x =seq_len(100), test = get_mean3(sim, phi1, phi2))
}
Graph the latter two series and compare them.
ARIMA11<-create_ARIMA_frame(sim_ARIMA11,.6,.6)
AR2<-create_AR2_frame(sim_ARIMA11,-.8,.3)
dat<-rbind(ARIMA11,AR2)
ggplot(dat, aes(x = x, y = test, color = model))+geom_line()
It seems that the AR model has far more fluctuations and the mean seems to be lower. This is confirmed by our prior knowledge that an AR with phi<0 will tend to oscillate around the mean.
9.7
Consider aus_airpassengers, the total number of passengers (in millions) from Australian air carriers for the period 1970-2011.
a. Use ARIMA() to find an appropriate ARIMA model. What model was selected. Check that the residuals look like white noise. Plot forecasts for the next 10 periods.
First lets look at our data.
air<-aus_airpassengers%>%
filter(Year<=2011)
air%>%autoplot()
## Plot variable not specified, automatically selected `.vars = Passengers`
Our data looks like it is non-seasonal and has a trend. Lets follow the following modeling procedure
Procedure:
1. Plot and identify any unusual Observations.
2. Transform Data if Necessary
3. if data are non-stationary, take first differences of data until stationary
4. Examine ACF/PACF
-is ARIMA(p,d,0) pr ARIMA(0,d,q) more appropriate?
5. Try chosen model(s) and use AIC_c to search for a better model
6. Check residuals from chosen model by plotting ACF of residuals, and doing a portmanteau test. If they do not look like white noise, try a modified model.
Plot and identify any unusual Observations.
After plotting above, there do not seem to be any outliers of note.
Transform Data if Necessary
Lets run a box cox and look at the suggested lambda
lambda<-BoxCox.lambda(air$Passengers)
lambda
## [1] 0.8927667
air%>%gg_tsdisplay((Passengers), plot_type = 'partial')
This is fairly close to 1, which would reccomend no box cox. Lets try modeling off of our untransformed data. Our initial ACF and PACF plots (undifferenced) also don’t make us think we need to transform the data.
if data are non-stationary, take first differences of data until stationary
The data clearly has trend, so lets use a unitroot test to confirm that we should take a difference.
air%>%
features(Passengers, unitroot_ndiffs)
## # A tibble: 1 x 1
## ndiffs
## <int>
## 1 2
This does confirm that a first difference is required.
air_t<-air%>%mutate(Passengers= difference(Passengers))
print(air_t%>%features(Passengers, unitroot_ndiffs))
## # A tibble: 1 x 1
## ndiffs
## <int>
## 1 1
air_t%>%gg_tsdisplay(Passengers, plot_type = 'partial')
## Warning: Removed 1 row(s) containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_point).
Looking at the variance in the difference, I am thinking maybe we do want to transform the data, lets check what that plot looks like
transformed_air_t<-air%>%
mutate(Passengers = forecast::BoxCox(Passengers, .5))
transformed_air_t%>%gg_tsdisplay(difference(Passengers), plot_type = 'partial')
## Warning: Removed 1 row(s) containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_point).
After transforming, the variance looks a lot more consistent, so we’ll use this data moving forward.
Examine ACF/PACF ,is ARIMA(p,d,0) or ARIMA(0,d,q) model appropriate?
There are no spikes outside of our critical values, so I assume that neither ARIMA is suggested?
Our PACF for the undifferenced data had lag 1 above our critical value, so I will include an AR(1) (ARIMA(1,1,0)) out of interest.
- Try chosen model(s) and use AIC_c to search for a better model
fit<-transformed_air_t%>%
model(
arima110 = ARIMA(Passengers~pdq(1,1,0)),
arima010 = ARIMA(Passengers~pdq(0,1,0)),
arima010_no_drift = ARIMA(Passengers~0+pdq(0,1,0), ),
arima212 = ARIMA(Passengers~pdq(2,1,2)),
arima212_no_c = ARIMA(Passengers~0+pdq(2,1,2)),
arima021 = ARIMA(Passengers~pdq(0,2,1)),
arima021_no_c = ARIMA(Passengers~0+pdq(2,1,2)),
search = ARIMA(Passengers)
)
## Warning: 1 error encountered for arima212_no_c
## [1] non-stationary AR part from CSS
## Warning: 1 error encountered for arima021_no_c
## [1] non-stationary AR part from CSS
a<-fit%>%pivot_longer(everything(), names_to = 'model_name', values_to = 'orders')
print(a$orders)
## <lst_mdl[8]>
## [1] <ARIMA(1,1,0) w/ drift> <ARIMA(0,1,0) w/ drift> <ARIMA(0,1,0)>
## [4] <ARIMA(2,1,2) w/ drift> <NULL model> <ARIMA(0,2,1)>
## [7] <NULL model> <ARIMA(0,1,0) w/ drift>
glance(fit)%>%arrange(AICc)%>%select(.model:BIC)
## # A tibble: 6 x 6
## .model sigma2 log_lik AIC AICc BIC
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 arima010 0.133 -16.3 36.6 36.9 40.0
## 2 search 0.133 -16.3 36.6 36.9 40.0
## 3 arima110 0.135 -16.2 38.3 39.0 43.5
## 4 arima021 0.141 -18.1 40.3 40.6 43.6
## 5 arima212 0.140 -15.3 42.7 45.1 52.9
## 6 arima010_no_drift 0.195 -24.7 51.4 51.5 53.1
Check residuals from chosen model by plotting ACF of residuals, and doing a portmanteau test. If they do not look like white noise, try a modified model.
fit%>%
select(search)%>%
gg_tsresiduals()
the residuals look like white noise, so we will select the random walk model.
b. Write the model in terms of the backshift operator.
ARIMA (0,1,0) in backshift notation:
$$
By_t = y_{t-1}\ and\ ARIMA(0,1,0) = y_t- y_{t-1}\ so\ ARIMA(0,1,0) = y_t-By_t = (1-B)y_t
$$
c. Plot forecasts from an ARIMA(0,1,0) model with drift and compare these to part a.
fit%>%
forecast(h=5)%>%
filter(.model =='arima010')%>%
autoplot(transformed_air_t)+
labs(title = 'ARIMA(0,1,0) Forecast')
This was the model I selected, so the forecasts are the same.
d. Plot forecasts from an ARIMA(2,1,2) model with drift and compare these to parts a and c. Remove the constant and see what happens.
fit%>%
forecast(h=5)%>%
filter(.model =='arima212')%>%
autoplot(transformed_air_t)+
labs(title = 'ARIMA(2,1,2) Forecast')
The removal of the constant results in a null-model. So you cannot plot the forecast. e. Plot forecasts from an ARIMA(0,2,1) model with a constant. What happens?
fit%>%
forecast(h=5)%>%
filter(.model =='arima021')%>%
autoplot(transformed_air_t)+
labs(title = 'ARIMA(0,2,1) Forecast')
When modeling without a constant, I get a null-model, but with a constant I get a prediction. ## 9.8
For the United States GDP series (from global_economy):
gdp<-global_economy%>%
filter(Country=='Germany')%>%
select(Year, GDP)
a. if necessary, find a suitable Box-Cox transformation for the data;
Examining data and testing the suggested lambda
gdp%>%gg_tsdisplay(difference(GDP), plot_type = 'partial')
## Warning: Removed 11 row(s) containing missing values (geom_path).
## Warning: Removed 11 rows containing missing values (geom_point).
It does look like the variance is not consistent, so I will apply the appropriate transformation.
lambda<-forecast::BoxCox.lambda(gdp$GDP)
gdp<-gdp%>%
mutate(GDP = BoxCox(GDP, lambda))
gdp%>%gg_tsdisplay(difference(GDP), plot_type = 'partial')
## Warning: Removed 11 row(s) containing missing values (geom_path).
## Warning: Removed 11 rows containing missing values (geom_point).
b. fit a suitable ARIMA model to the transformed data using ARIMA();
Judging by the ACF and PACF plots, I wouldn’t reccomend any particular model. The values above the critical values are too far out to use as values for p or q. I will include ARIMA(0,1,3) based on the untransformed plot just to check.
c. try some other plausible models by experimenting with the orders chosen;
fit<-gdp%>%
model(
arima110 = ARIMA(GDP~pdq(1,1,0)),
arima010 = ARIMA(GDP~pdq(0,1,0)),
arima013 = ARIMA(GDP~pdq(0,1,3)),
arima212 = ARIMA(GDP~pdq(2,1,2)),
arima021 = ARIMA(GDP~pdq(0,2,1)),
search = ARIMA(GDP)
)
a<-fit%>%pivot_longer(everything(), names_to = 'model_name', values_to = 'orders')
print(a$orders)
## <lst_mdl[6]>
## [1] <ARIMA(1,1,0) w/ drift> <ARIMA(0,1,0) w/ drift> <ARIMA(0,1,3) w/ drift>
## [4] <ARIMA(2,1,2) w/ drift> <ARIMA(0,2,1)> <ARIMA(0,1,1) w/ drift>
glance(fit)%>%arrange(AICc)%>%select(.model:BIC)
## # A tibble: 6 x 6
## .model sigma2 log_lik AIC AICc BIC
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 search 17661. -300. 607. 607. 613.
## 2 arima021 19908. -302. 607. 607. 611.
## 3 arima110 18165. -301. 608. 608. 614.
## 4 arima010 18938. -302. 609. 609. 613.
## 5 arima013 18110. -300. 610. 611. 620.
## 6 arima212 17518. -299. 610. 612. 622.
a$orders[6]
## <lst_mdl[1]>
## [1] <ARIMA(0,1,1) w/ drift>
The searched model is ARIMA(0,1,1) w/drift , or simple exponential smoothing with growth.
d. choose what you think is the best model and check the residual diagnostics;
fit%>%
select(search)%>%
gg_tsresiduals()
## Warning: Removed 10 row(s) containing missing values (geom_path).
## Warning: Removed 10 rows containing missing values (geom_point).
## Warning: Removed 10 rows containing non-finite values (stat_bin).
The residuals look normal without trend in the innovation residuals.
e. produce forecasts of your fitted model. Do the forecasts look reasonable?
fit%>%
forecast(h=5)%>%
filter(.model =='search')%>%
autoplot(gdp)+
labs(title = 'ARIMA(0,2,1) Forecast')
## Warning: Removed 10 row(s) containing missing values (geom_path).
This looks like a fine prediction. It generally follows the trend of the data. f. compare the results with what you would obtain using ETS() (with no transformation).
global_economy %>%
filter(Country == 'Germany') %>%
select(Year, GDP)%>%
na.omit()%>%
model(
ETS(GDP)
) %>%
forecast(h = 5) %>%
autoplot(global_economy %>%
filter(Country == 'Germany') %>%
select(Year, GDP))
## Warning: Removed 10 row(s) containing missing values (geom_path).
Compared to the ETS model, The main difference is the spread of the confidence intervals. It appears our ARIMA model gives similar point forecasts, but a higher degree of certainty.