9.1a Response

The autocorrelation plots differ because of 1. sampling variation in how acfs are calculated and 2. sample size differences.

Because white noise is comprised of i.i.d observations their autocorrelations are subject to the law of large numbers. The autocorrelations should not exceed the confidence bands 95% of the time according to Section 2.9 on White Noise.

All 3 time series look like white noise because autocorrelations stay inside the band thresholds of the plots.

9.1b Response

Section 2.9 on White Noise states that 95% of the spikes in the ACF should lie within \(\pm 2/ \sqrt{T}\) where \(T\) is the length of the time series. Since we have three time series of lengths \(36, 360, 1000\), the size of the bands, i.e. \(T\), is different.

## [1] "T=36    band upper bound: 0.333"

## [1] "T=360   band upper bound: 0.105"

## [1] "T=1000  band upper bound: 0.063"

As the time series length increases, the confidence bands shrink.

9.3a Response

We observe that Turkey’s nominal GDP displays an exponential growth trend punctuated by business cycle recessions and booms.

We observe that the 1960-1961 change in log GDP is -55%. This means the economy completely collapsed. We regard this event as an outlier and remove it from the time series. We determined from Erik Zurcher’s Turkey: A Modern History, 4th edition (Chapters 13, 14):

• In 1960, Turkey experienced a military coup and constitutional crisis.
• Also, Turkey established its Second Turkish Republic which lasted from 1960-1980.
• Turkey accepted the first IMF bailout package which required it to liberalize trade policies, devalue the currency and adopt more pro-capitalist practices.

For these reasons, I choose to drop 1960 GDP from consideration and begin the series at 1961. Then the first differences are observed starting in 1962.

Next, we have to test for stationarity of the first differenced time series. We run two unit root tests: kpss and adf.

## # A tibble: 1 x 3
##   Country kpss_stat kpss_pvalue
##   <fct>       <dbl>       <dbl>
## 1 Turkey      0.100         0.1

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

## Warning in tseries::adf.test(turkey_gdp2$dlogGDP[2:57]): p-value smaller than
## printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  turkey_gdp2$dlogGDP[2:57]
## Dickey-Fuller = -4.5696, Lag order = 3, p-value = 0.01
## alternative hypothesis: stationary

We conclude based on the tests that \(\Delta(log(GDP))\) is likely stationarity as an approximation.

9.3b Response

We clearly observe seasonality, trend and non-constant variance in the above graph. The series appear to increase in volatility with time. Therefore, we will use a Box-Cox transformation to stabilize the variance before differencing to remove seasonality and trend. Using the Guerrero method we get a Box-Cox \(\lambda = -0.488\).

## [1] -0.04884781

Taking a further look with a gg_season plot, we see a lag 4 seasonality in the transformed series \(Y_t\) below:

To remove seasonality, we apply \(1-L^4\). To remove the trend, we apply \(1-L\). The resulting series \(X_t\) is defined below in terms of \(Y_t\)

\[X_t = (1-L^4)(1-L)Y_t = ( Y_{t} - Y_{t-4} ) - ( Y_{t-1} - Y_{t-5})\]

Plotting it, we see that \(X_t\) looks stationary. Does it pass the statistical tests?

The KPSS test shows \(X_t\) is stationary:

## # A tibble: 1 x 3
##   State    kpss_stat kpss_pvalue
##   <chr>        <dbl>       <dbl>
## 1 Tasmania    0.0467         0.1

The ADF test shows \(X_t\) to be stationary with k=6 lags. (By default \(k\) is lower, so we are running a stricter test.)

## Warning in tseries::adf.test(tas_deseasoned$X[!is.na(tas_deseasoned$X)], : p-
## value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  tas_deseasoned$X[!is.na(tas_deseasoned$X)]
## Dickey-Fuller = -4.1507, Lag order = 6, p-value = 0.01
## alternative hypothesis: stationary

9.3c Response

We observe in the raw data plot of souvenirs sales a pattern with:

• Annual spikes in December (which is summer in Australia). So a distinct pattern of seasonality
• Increasing trend
• Increasing variance

It appears that the detrended, deseasoned log Sales time series is stationary according to the KPSS test. Specifically, if \(X_t\) is the time series that we test and \(Y_t\) is the Sales:

\[X_t = (1-L)(1-L^{12})Y_t\]

## # A tibble: 1 x 2
##   kpss_stat kpss_pvalue
##       <dbl>       <dbl>
## 1    0.0381         0.1

We confirm that by running the ADF test with \(k=6\) lags:

## 
##  Augmented Dickey-Fuller Test
## 
## data:  souv$dS[!is.na(souv$dS)]
## Dickey-Fuller = -4.7674, Lag order = 6, p-value = 0.01
## alternative hypothesis: stationary

9.6a Response

We run the code as provided to simulate the AR(1) model.

set.seed(100)

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)

The simulated output is stored is a tsibble called sim.

9.6b Response

Plotting the series we get:

We experiment with plotting AR(1) time series using four values of \(\phi_1=0.05, 0.6, 0.99, 1.05\).

The same white noise process \(\epsilon_t\) is used for all four plots below.

set.seed(100)

y <- numeric(100)
w <- numeric(100)
u <- numeric(100)
z <- numeric(100)

e <- rnorm(100)

for(i in 2:100)
{
  y[i] <- 0.6 *  y[i-1] + e[i]
  w[i] <- 0.05 * w[i-1] + e[i]
  u[i] <- 0.99 * u[i-1] + e[i]
  z[i] <- 1.05 * z[i-1] + e[i]
}
sim <- tsibble(idx = seq_len(100), y = y, w = w, u = u, z = z, index = idx)

We observe that:

• \(\phi_1=0.05\) produces a white noise like time series. The small contribution of the prior series value \(X_{t-1}\) means the noise term is dominant.
• \(\phi_1=0.60\) produces a series that has some autoregressive properties.
• \(\phi_1=0.99\) produces a series which is almost like a random walk. Since the contribution of the prior series value \(X_{t-1}\) is .99, the past history has a very significant impact on the current level.
• \(\phi_1=1.05\) causes the series to exponentially explode.

9.6c Response

We are asked to generate an MA(1) model with \(\phi_1=0.5\) and \(\sigma^2=1\).

set.seed(100)
y <- numeric(100)
e <- rnorm(100)
for(i in 2:100)
  y[i] <- 0.6*e[i-1] + e[i]
sim <- tsibble(idx = seq_len(100), y = y, index = idx)

9.6d Response

Let’s try \(\theta_1 = -1, -\frac{1}{2}, 0.05, \frac{1}{2}, 0.99, 1.3\) for the \(MA(1)\) model:

set.seed(100)

y <- numeric(100)
w <- numeric(100)
u <- numeric(100)
z <- numeric(100)
r <- numeric(100)
s <- numeric(100)

e <- rnorm(100)

thetas = c( -1.0, -0.5, 0.05 ,  0.5, 0.99 , 1.3 )

for(i in 2:100)
{
  u[i] <- thetas[1] * e[i-1] + e[i]
  w[i] <- thetas[2] * e[i-1] + e[i]
  y[i] <- thetas[3] * e[i-1] + e[i]
  z[i] <- thetas[4] * e[i-1] + e[i]
  r[i] <- thetas[5] * e[i-1] + e[i]
  s[i] <- thetas[6] * e[i-1] + e[i]
}
sim <- tsibble(idx = seq_len(100),   u = u, w = w, y = y, z = z, r = r, s = s,  index = idx)

##        u        w        y        z        r        s 
## 1.618886 1.254229 1.007038 1.009586 1.224021 1.433917

We observe that:

• \(\theta_1 = -1\) causes extreme oscillations between consecutive observations. The entire series is very choppy.
• \(\theta_1 = -0.5\) also shows extreme oscillations between consecutive observations, but the series shows less extremem output.
• \(\theta_1 = 0.05\) shows damped amplitudes.
• \(\theta_1 = 0.5\) shows moderate amplitudes and smoother oscillations between consecutive observations.
• \(\theta_1 = 0.99, 1.3\) show larger amplitudes but even smoother oscillations. Locally sustained trends, peaks and valleys become more visible.

9.6e Response

The plot below shows an \(ARMA(1,1)\) model realization with \(\theta_1 = 0.6\) and \(\phi_1 = 0.6\).

\[Y_t = 0.6 \cdot Y_{t-1} + 0.6 \cdot \epsilon_{t-1} + \epsilon_{t}\]

set.seed(100)
y <- numeric(100)
e <- rnorm(100)

for(i in 2:100)
  y[i] <- 0.6 * y[i-1] + 0.6*e[i-1] + e[i]
sim <- tsibble(idx = seq_len(100), y = y, index = idx)

ARMA(1,1) realization

9.6f Response

We generate an AR(2) model with non-stationary coefficients.

\[X_{t} = -0.8 \cdot X_{t-1} + 0.3 \cdot X_{t-2} + \epsilon_{t}\]

set.seed(100)
zz <- numeric(100)
ez <- rnorm(100)

for(i in 3:100)
  zz[i] <- -0.8 * zz[i-1] + 0.3*zz[i-2] + ez[i]
simz <- tsibble(idx = seq_len(100), zz = zz, index = idx)

AR(2) realization

9.6g Response

The AR(2) model has non-constant variance, no trend (with zero mean). Its growth appears to be exponential because the amplitudes of the oscillations grow very fast.

In contrast, the ARMA(1,1) model has constant variance and appears to have no trend. It exists some trends consists with autocorrelation properties.

9.7a Response

Before using the ARIMA() function on the raw data series, we do exploratory data analysis by visualizing the air passenger time series in 4 plots labeled P1-P4.

⊕The light color plots P1, P3 use the raw data series: Aus Air Passengers in Millions from 1970-2016. P1 shows raw data overlaid by a red trend line using the STL decomposition method. P3 shows the STL decomposition components. The dark themed plots P2, P4 use the log transformed data series: log(Aus Air Passengers in Millions). Their purpose is similar to P1 and P3 respectively.

Observations:

• In P1, the raw data and trend line are clearly increasing but possibly polynomial or exponential.
  
• In P3, the residuals are increasing in amplitude as time passes. This suggests non-constant variance.
• In P2, the log data series clearly suggest a linear trend. This implies an exponential growth trend of the raw time series. By contrast, if P1 has displayed polynomial growth. then P2 would display a concave trend that increases but bends down.
• In P4, the residuals show stable variance.

These observations suggest a good ARIMA model should forecast exponential future growth.

Returning to the task, the ARIMA() function selects an ARIMA(0,2,1) model as its optimal choice below.

fit_default = ausp %>% model(ARIMA(Passengers))
report(fit_default)

## Series: Passengers 
## Model: ARIMA(0,2,1) 
## 
## Coefficients:
##           ma1
##       -0.8963
## s.e.   0.0594
## 
## sigma^2 estimated as 4.308:  log likelihood=-97.02
## AIC=198.04   AICc=198.32   BIC=201.65

The KPSS test suggests that the twice-differenced Passengers time series is stationary.

## # A tibble: 1 x 2
##   kpss_stat kpss_pvalue
##       <dbl>       <dbl>
## 1    0.0361         0.1

We plot the diagnostics first.

• ACF plot shows a white noise like pattern of autocorrelations.
• the histogram of innovation residuals shows some outliers on the right tail due to the 2010 spike.
  
• Residuals volatility seems to grow over time

The innovation residuals pass the Ljung Box test. We conclude they look like white noise.

## # A tibble: 1 x 3
##   .model            lb_stat lb_pvalue
##   <chr>               <dbl>     <dbl>
## 1 ARIMA(Passengers)    6.70     0.754

Lastly, the 10 year model forecasts are shown below. We observe that the mean growth is an increasing straight line. It appears to simply extrapolate the tangent line defined by the derivative at the end point. Thus, the default ARIMA model produces linear not exponential growth forecasts.

9.7b Response

The above “default” fitted ARIMA model ARIMA(0,2,1) can be written in terms of the lag (a.k.a. backshift) operator. Let \(Y_t\) denote the raw passengers time series.

The left hand side shows the degree of differencing: 2. There is no AR(p) component otherwise the AR(p) terms goes on the LHS.

The right hand side shows the order of the MA(q) model: 1. So we get:

\[ (1-B)^2 Y_t = (1 + \theta_1B)\epsilon_t \] This can be rewritten to isolate \(Y_t\) on the LHS:

\[ \begin{eqnarray} (1-2B +B^2)Y_t &=& \epsilon_t + \theta_1 B\epsilon_t \\ Y_t -2Y_{t-1}+Y_{t-2} &=& \epsilon_t + \theta_1 \epsilon_{t-1}\\ Y_t &=& 2Y_{t-1} - Y_{t-2} + \theta_{1}\epsilon_{t-1} + \epsilon_{t} \\ Y_t &=& 2Y_{t-1} - Y_{t-2} - 0.8963 \epsilon_{t-1} + \epsilon_{t} \\ \end{eqnarray} \]

9.7c Response

We will fit an ARIMA(0,1,0) model with drift below. We use the log(Passengers) as the input variable. This model is basically equivalent to an exponential model on Passengers estimated at discrete time steps.

fit_log = ausp %>% model(ARIMA(log(Passengers) ~ 1 + pdq(0,1,0) ) )
report(fit_log)

## Series: Passengers 
## Model: ARIMA(0,1,0) w/ drift 
## Transformation: log(Passengers) 
## 
## Coefficients:
##       constant
##         0.0499
## s.e.    0.0104
## 
## sigma^2 estimated as 0.005089:  log likelihood=56.69
## AIC=-109.38   AICc=-109.1   BIC=-105.72

We see that the transformed series passes the unit root test for stationarity below.

ausp %>% mutate(diff_logpass = difference(log(Passengers) )) %>%
   features( diff_logpass, unitroot_kpss )

## # A tibble: 1 x 2
##   kpss_stat kpss_pvalue
##       <dbl>       <dbl>
## 1    0.0516         0.1

We also see that the residual diagnostics look good.

And the model residuals pass the white noise test.

## # A tibble: 1 x 3
##   .model                                    lb_stat lb_pvalue
##   <chr>                                       <dbl>     <dbl>
## 1 ARIMA(log(Passengers) ~ 1 + pdq(0, 1, 0))    8.19     0.610

Finally, plotting the forecast of this model produces the following plot:

We draw comparisons with the previous model in a.

• Whereas the ARIMA(0,2,1) shows linear growth, the ARIMA(0,1,0) model shows the desired exponential growth.
• Both models satisfy the unit root and white noise requires.
• Whereas the ARIMA(0,2,1) model residuals show increasing variance (heteroskedasticity). The ARIMA(0,1,0) model residuals exhibit the desired homoskedasticity.

9.7d Response

Now we consider an ARIMA(2,1,2) model with drift on the log(Passengers).

fit_arima212c = ausp %>% model(ARIMA(log(Passengers) ~ 1 + pdq(2,1,2) ) )
report(fit_arima212c)

## Series: Passengers 
## Model: ARIMA(2,1,2) w/ drift 
## Transformation: log(Passengers) 
## 
## Coefficients:
##          ar1      ar2      ma1     ma2  constant
##       0.8362  -0.2286  -1.0767  0.0767    0.0197
## s.e.  0.8981   0.5910   0.9354  0.9333    0.0006
## 
## sigma^2 estimated as 0.004391:  log likelihood=60.83
## AIC=-109.67   AICc=-107.52   BIC=-98.7

Since we already know the \(\Delta(log(Passengers))\) is stationary from part c, we just need to verify the residuals look like white noise with the Ljung Box test. They pass that test.

## # A tibble: 1 x 3
##   .model                                    lb_stat lb_pvalue
##   <chr>                                       <dbl>     <dbl>
## 1 ARIMA(log(Passengers) ~ 1 + pdq(2, 1, 2))    6.70     0.754

The visual inspection of diagnostic plots below concurs with the Ljung-Box test result.

Now, we plot the resulting forecasts below.

Let’s compare the model outputs in d. with the models from a. and c.

• the ARIMA(2,1,2) and ARIMA(0,1,0) models in d. and c. produce similar forecasts as shown below.
• their model residuals pass the white noise test and the input series passes unit root tests.
• their forecasts appear non-linear - likely exponential – and are better than the model in a.

If we remove the constant drift term and fit the ARIMA(2,1,2) model, we get an error message that the resulting model has a non-stationary AR component. So the resulting model cannot be used for forecasting.

fit_arima212 = ausp %>% model(ARIMA(log(Passengers) ~ 0 + pdq(2,1,2) ) )

## Warning: 1 error encountered for ARIMA(log(Passengers) ~ 0 + pdq(2, 1, 2))
## [1] non-stationary AR part from CSS

report(fit_arima212)

## Series: Passengers 
## Model: NULL model 
## Transformation: log(Passengers) 
## NULL model

9.7e Response

We observe that ARIMA(0,2,1) model introduces a quadratic or higher order polynomial trend. This is discouraged by the model authors.

fit_arima021c = ausp %>% model(ARIMA(Passengers ~ 1 + pdq(0, 2, 1) ) )

## Warning: Model specification induces a quadratic or higher order polynomial trend. 
## This is generally discouraged, consider removing the constant or reducing the number of differences.

report(fit_arima021c)

## Series: Passengers 
## Model: ARIMA(0,2,1) w/ poly 
## 
## Coefficients:
##           ma1  constant
##       -1.0000    0.0571
## s.e.   0.0585    0.0213
## 
## sigma^2 estimated as 3.855:  log likelihood=-95.1
## AIC=196.21   AICc=196.79   BIC=201.63

If we look at the 10 year ahead forecast, the mean forecast is an upward trending non-linear pattern.

However, the growth rate is slower than that of our models in parts c and d. While this model produces quadratic growth, models in parts c. and d. produce exponential growth.

9.8a Response

A suitable Box-Cox transformation is to take the log of GDP, i.e. use \(\lambda=0\). However, we investigate the use of the guerrero algorithm to find an optimal Box-Cox parameter. It produces \(\lambda=0.28\).

• When we plot below the two Box-Cox transformation candidates side by side, there is no discernible benefit of \(\lambda=0.28\) in stabilizing variance of GDP over \(\lambda=0\).
  
• The benefit of using the \(\lambda=0\) option is that \(\Delta(log(GDP))\) approximates the annual GDP growth rate, an extremely important economic indicator.
  
• The Box-Cox transform hides the slowdown in GDP growth in the last 25 years. This is shown as a slight concavity in the log(GDP) plot but the bend in the curve is straighter in the Box-Cox 0.28 plot.

## [1] 0.2819443

Therefore, we will simply use the log transformation in the rest of this exercise.

9.8b

To fit a suitable ARIMA model, we first perform some EDA, test for stationarity and use a search process to fit an ARIMA(1,1,0) model on log(GDP).

⊕Plots P1, P3, P5 show US GDP in raw form, in log form and annual changes of log. Plots P2, P4, P6 show the equivalent output for the US GDP per capita. GDP per capita is obtained by dividing nominal GDP by Population.

Based on the EDA above, we conclude:

• \(\Delta(log(GDP))\) (Plot P5) or \(\Delta(log(GDP per capita))\) (Plot P6) could be stationary series. They show little or no significant trend.

• US GDP, US GDP per capita and their logs are both choices. Plots P1-P4 are ruled out.

Next, we will see the Unit Root tests fail for US GDP over the long term, i.e. 1960-2017. This problem is well known in the econometric literature. The failures of the ADF and KPSS tests on the full time series are shown below.

In fact, modern GDP forecasting is done by the Federal Reserve bank of New York and Atlantic using a completely different approach. They use a state space model approach using a dynamic factor model that takes around 30 predictor variables. These are received in real-time as news releases which update the GDP forecast daily. The model expects information to arrive on a sporadic lagging basis. The model builders objective is to produce 1 quarter to 2 year ahead forecasts. In contrast, our dataset and forecast horizon is 10 years or longer.

Links to their model and approach are here: (https://www.newyorkfed.org/research/policy/nowcast/faqs.html)

## # A tibble: 1 x 3
##   Country       kpss_stat kpss_pvalue
##   <fct>             <dbl>       <dbl>
## 1 United States     0.963        0.01

## 
##  Augmented Dickey-Fuller Test
## 
## data:  usgdp$dlogGDP[2:length(usgdp$dlogGDP)]
## Dickey-Fuller = -2.6551, Lag order = 12, p-value = 0.3103
## alternative hypothesis: stationary

A solution is to truncate the data series to observations after 1995. We can see that the KPSS test confirms stationarity of the truncated series. I would add that the ADF test FAILS. So this discrepancy in unit root tests will require deeper consideration out of the scope of this assignment. For now, I model the series as if \(\Delta(log(GDP_t))\) is stationary.

usgdp2 = usgdp %>% filter(Year > 1995)

usgdp2 %>% features( dlogGDP, unitroot_kpss)

## # A tibble: 1 x 3
##   Country       kpss_stat kpss_pvalue
##   <fct>             <dbl>       <dbl>
## 1 United States     0.334         0.1

Now we fit an ARIMA model with a drift and constrain \(d=1\) the order of differencing on the \(X_t = log(GDP)\). We use a search process to allow \(p=0 ... 12\) and \(q=0..5\).

fit_usgdp = usgdp2 %>% model( m1 = ARIMA( log(GDP) ~ 1 + pdq(p = 0:12, d=1, q = 0:5 ) ) )

report(fit_usgdp)

## Series: GDP 
## Model: ARIMA(1,1,0) w/ drift 
## Transformation: log(GDP) 
## 
## Coefficients:
##          ar1  constant
##       0.4751    0.0222
## s.e.  0.1885    0.0035
## 
## sigma^2 estimated as 0.0003544:  log likelihood=56.02
## AIC=-106.04   AICc=-104.63   BIC=-102.91

ARIMA produces an ARIMA(1,1,0) model on \(X_t\) whose equation is:

\[\begin{eqnarray} (1 - \phi_1 B)(1 - B)X_t&=& c + \epsilon_t\\ (1 - \phi_1 B - B + \phi_1B^2) X_t & = & c + \epsilon_t \\ X_t& = & c + \epsilon_t + (1 + \phi_1)X_{t-1} - \phi_1 X_{t-2}\\ X_t& = & 0.0222 + \epsilon_t + (1 + .4751)X_{t-1} - 0.4751X_{t-2} \\ X_t& = & 0.0222 + 1.4751X_{t-1} - 0.4751 X_{t-2} + \epsilon_{t} \end{eqnarray} \]

9.8c Response

To find other models, we look at the ACF and PACF of the series \(Y_t = \Delta X_t = log(GDP_t) - log(GDP_{t-1})\). The series diagnostic plot suggests that acf and pacf has significant values at lags 1 only.

Consequently, we consider 3 candidate models based on ACF and PACF information.

We specify them as ARIMA models on the same variable \(X_t\). We include a drift term because the US economy has a positive growth rate with deviations during economic downturns or recessions.

• Model m1 (the preferred one) is stated below in part b, but repeated here. It is ARIMA(1, 1, 0) with drift.
• Model m2 is ARIMA(1, 1, 1) with drift.
• Model m3 is ARIMA(2, 1, 2) with drift.

fit_usgdp = usgdp2 %>% model( m1 = ARIMA( log(GDP) ~ 1 + pdq(p = 0:12, d=1, q = 0:5 ) ) ,
                              m2 = ARIMA( log(GDP) ~ 1 + pdq(p = 1, d = 1, q = 1 ) ) ,
                              m3 = ARIMA( log(GDP) ~ 1 + pdq(p = 2, d = 1, q = 2 ) )
                              )

Their parameters are reported below for completeness.

Model m1 has the following parameters (repeating the output from part b.)

## Series: GDP 
## Model: ARIMA(1,1,0) w/ drift 
## Transformation: log(GDP) 
## 
## Coefficients:
##          ar1  constant
##       0.4751    0.0222
## s.e.  0.1885    0.0035
## 
## sigma^2 estimated as 0.0003544:  log likelihood=56.02
## AIC=-106.04   AICc=-104.63   BIC=-102.91

Model m2 has the following parameters. The coefficients values are plausible but the AICc is inferior.

## Series: GDP 
## Model: ARIMA(1,1,1) w/ drift 
## Transformation: log(GDP) 
## 
## Coefficients:
##          ar1     ma1  constant
##       0.2884  0.2600    0.0300
## s.e.  0.3499  0.3349    0.0044
## 
## sigma^2 estimated as 0.0003652:  log likelihood=56.29
## AIC=-104.58   AICc=-102.08   BIC=-100.4

Model m3 has the following parameters. The AR and the MA coefficients are unstable because they are opposite in sign and same magnitude. For example, \(\phi_1 = 1.11\) while \(\phi_2 = -0.948\).

## Series: GDP 
## Model: ARIMA(2,1,2) w/ drift 
## Transformation: log(GDP) 
## 
## Coefficients:
##          ar1      ar2      ma1     ma2  constant
##       1.1101  -0.9487  -0.9414  0.9999    0.0346
## s.e.  0.2065   0.1009   0.3803  0.5288    0.0034
## 
## sigma^2 estimated as 0.0003386:  log likelihood=56.92
## AIC=-101.85   AICc=-95.85   BIC=-95.58

9.8d Response

We choose that model m1 (the ARIMA(1,1,0)) as our best candidate because it:

• is parsimonious.
• gives similar predictions as the more complex models.
• has the lowest AICc score

We also check their residual diagnostics in the following charts. None of the model residuals look like white noise, but there are only 23 data points. However, the acf graphs don’t significant significant autocorrelations.

9.8e Response

We see that all 3 ARIMA models produce similar forecasts but different confidence intervals. All 3 forecasts look reasonable.

In general, model m3 produces the narrowest confidence intervals. Model m2 produced the widest intervals. Model m1, our preferred one, produces intermediate confidence intervals.

The forecasts are around 4.3% annualized growth which is somewhat plausible as the current GDP growth forecast by the New York Fed Nowcast model is 3.8% for Q3 2021. (https://www.newyorkfed.org/research/policy/nowcast.html) However, in the longer term, GDP forecast show be in the range of 1-4 percent with a 2-3 percent midrange.

US GDP Forecast Trillions

9.8f Response

## Series: GDP 
## Model: ETS(M,A,N) 
##   Smoothing parameters:
##     alpha = 0.9998997 
##     beta  = 0.0001002624 
## 
##   Initial states:
##          l[0]         b[0]
##  7.545823e+12 527283884848
## 
##   sigma^2:  3e-04
## 
##      AIC     AICc      BIC 
## 1225.410 1229.160 1230.865

In the above chart showing 10 year ahead forecasts from the ARIMA(1,1,0) model (called m1) versus the exponential smoothing model ETS(M,A,N), we see that

• The ETS model produces linear GDP growth
• The ARIMA(1,1,0) model predicts much higher growth with an exponential trend.

We do not find the ETS model useful in forecasting GDP.