Model identification

Use the following command to load two data sets into your R session as the ts-objects named gnp and birth.

load(url("https://www.math.ntnu.no/emner/TMA4285/2017h/project2.RData"))

The data sets contain quarterly data on the U.S. grand national product between 1947 and 2002 and monthly birth data in the US from 1948 to 1978.

  1. Select and fit an ARIMA model to the GNP data. This includes possible transformations of the data and seasonal models if needed. Compute a 10-year forecast.

  2. Select and fit an ARIMA model for the birth data. Again, you may transform the data and use seasonal models. Compute a 12-month ahead forecast.

You will need to use functions acf, pacf, diff, tsdiag, arima, and predict.Arima in R.

State space models and the Kalman filter

Suppose that \(\{X_t\}\) follows a second order autoregressive process, that is, \[ (1-\phi_1 B -\phi_2 B^2) X_t = a_t \tag{1} \] but that we only observe \(X_t\) indirectly through the observation equation \[ Z_t = X_t + b_t \tag{2} \] where \(a_t\) and \(b_t\) are independent Gaussian white noise processes with variances \(\sigma_a^2\) and \(\sigma_b^2\).

  1. Show that \(Z_t\) follows an ARIMA\((2,0,2)\) process. Hint: Apply \(\phi(B)\) to both sides of (2). What can you say about the autocovariance function of the right hand side of the resulting equation?

Suppose that we the following observations of \(Z_1,Z_2,\dots,Z_n\), \(n={50}\). 

c(-7.59, -8.81, -6.13, -5.7, -4.23, 3.21, 8.94, 13.75, 19.03, 
17.37, 15.76, 8.96, 3.49, -2.9, -8.43, -13.21, -14.8, -18.18, 
-12.86, -7.72, -2.5, 3.86, 10.07, 12.26, 16.75, 16.37, 14.06, 
6.33, -0.11, -7.98, -12.51, -16.1, -15.53, -12.95, -11.32, -6.33, 
-2.15, 0.66, 4, 8.52, 11.79, 14.84, 14.09, 10.13, 9.82, 6.61, 
4.72, 2.58, 1.04, -0.37)

  1. Fit a ARIMA\((2,0,2)\) to the data and compute the \(l\)-step ahead forecast and associated forecast limits (using predict.Arima) for \(l=1,2,\dots,10\).

  2. Write (1) and (2) in state space form. Implement an R function kalman carrying out the Kalman forecasting and filtering recursion. Your function should take the model parameters \(\phi_1,\phi_2,\sigma_a^2,\sigma_b^2\) as a first vector argument and the data \(z_1,\dots,z_n\) as a second vector argument. What’s the joint distribution of \(X_1,X_2\)? (see lecture summary, pp. 8-10). Your function should return a list containing estimates of the state vector (the mean vector and its covariance matrix) after each forecast and filtering step at each time point \(t=1,2,\dots,n\). Suiteable data structures for storing these vectors and matrices in R are \(2\times n\) and \(2\times 2\times n\) arrays (e.g. array(dim=c(2,2,n)))).

  3. Run the filter based on what you think are reasonable parameter values and plot the results together with the observed data to check that your implementation of the algorithm is working.

  4. Modify your R function such that it also computes and optionally (if an addition argument nll=TRUE) instead returns the negative log likelihood only. Maximise the likelihood by minimising the negative log likelihood using by doing optim(startvals, kalman, z=z, nll=TRUE). You may need to play around with different starting values for the parameters. Compute approximate standard errors of the estimates by inverting (R function solve) the observed fisher information matrix (returned by optim if the optional argument hessian=TRUE).

  5. To check that your code works, rerun the filter based on the MLEs of the parameters and plot the filtered and forecasted states again. Briefly discuss if the results look reasonable.

  6. How do the estimates of \(\phi_1\) and \(\phi_2\) based on the the method in point 2 and 5 4 and 7 compare? In particular, is anything gained by working with the state space model rather than the ARIMA(2,0,2) model? Which is the most parsimonous model jugded by the AIC-values? (\(AIC = 2 p - 2 \log L\) where \(p\) is the number of estimated parameters and \(\log L\) is the maximum log likelihood under a given model).