Question 1

rm(list = ls(all = TRUE)) # Remove all objects in the memory
graphics.off()            # Close all graphs
library(tseries)
library(forecast)
library(astsa)
library(ggplot2)
work_dir = getwd()        # get work directory path
setwd(work_dir)           # set the work directory
column.names <- c("r", "money", "price", "y") # define the labels
usa.macro.monthly <- read.table("usa_macro_monthly.txt", header = FALSE, 
                                sep = "", col.names = column.names) # Load the data
attach(usa.macro.monthly)
gdp <- ts(y, start = c(1959, 1), end = c(1998, 12), frequency = 12)
gdp.sub <- window(gdp, start=c(1965, 1),end=c(1998, 12))
lngdp <- log(gdp.sub) # Log(GDP)
lngdp.d1 <- diff(lngdp, lag = 1) # first lag difference

(a)

autoplot(lngdp, main = "Log GDP")

autoplot(lngdp.d1, main = "First lag difference")

(b)

acf2(lngdp.d1, max.lag = 9)

##        ACF  PACF
##  [1,] 0.37  0.37
##  [2,] 0.27  0.15
##  [3,] 0.22  0.10
##  [4,] 0.14  0.01
##  [5,] 0.04 -0.06
##  [6,] 0.02 -0.02
##  [7,] 0.03  0.02
##  [8,] 0.06  0.06
##  [9,] 0.05  0.02

(c)

Estimation of ARMA Models

summary(arma.1 <- arima(lngdp.d1, order = c(1, 0, 1), method = "CSS"))
## 
## Call:
## arima(x = lngdp.d1, order = c(1, 0, 1), method = "CSS")
## 
## Coefficients:
##          ar1      ma1  intercept
##       0.7278  -0.4283     0.0025
## s.e.  0.0731   0.0953     0.0008
## 
## sigma^2 estimated as 5.28e-05:  part log likelihood = 1426.75
## 
## Training set error measures:
##                        ME        RMSE         MAE MPE MAPE      MASE
## Training set 6.047247e-06 0.007257602 0.005246083 NaN  Inf 0.7896541
##                      ACF1
## Training set -0.004430946
summary(arma.2 <- arima(lngdp.d1, order = c(2, 0, 2), method = "CSS"))
## 
## Call:
## arima(x = lngdp.d1, order = c(2, 0, 2), method = "CSS")
## 
## Coefficients:
##          ar1     ar2     ma1      ma2  intercept
##       0.1974  0.3888  0.1030  -0.2368     0.0024
## s.e.  0.4807  0.3512  0.4822   0.2200     0.0008
## 
## sigma^2 estimated as 5.278e-05:  part log likelihood = 1426.84
## 
## Training set error measures:
##                         ME        RMSE         MAE MPE MAPE      MASE
## Training set -2.779439e-06 0.007247068 0.005224524 NaN  Inf 0.7864091
##                      ACF1
## Training set -0.004427706

(d)

AIC and BIC

arma.1$loglik
## [1] 1426.754
arma.2$loglik
## [1] 1426.843

(e)

p1 <- c(1, -arma.2$coef[1:2]) # Obtain the characteristic equation of the model
p1 
##                   ar1        ar2 
##  1.0000000 -0.1974482 -0.3888066
Re(polyroot(p1))
## [1]  1.369798 -1.877629
p2 <- c(1, -arma.2$coef[3:4])
Re(polyroot(p2))
## [1] 0.2176059 0.2176059
plot.Arima(arma.2)

(f)

res.1 <- arma.1$residuals
res.2 <- arma.2$residuals

Ljung-Box test for serial correlation

Box.test(res.1, lag = 6, type = "Ljung-Box") 
Box.test(res.1, lag = 10, type = "Ljung-Box") 
Box.test(res.1, lag = 14, type = "Ljung-Box") 
Box.test(res.2, lag = 6, type = "Ljung-Box") 
Box.test(res.2, lag = 10, type = "Ljung-Box") 
Box.test(res.2, lag = 14, type = "Ljung-Box") 
## 
##  Box-Ljung test
## 
## data:  res.1
## X-squared = 4.0816, df = 6, p-value = 0.6656
## 
## 
##  Box-Ljung test
## 
## data:  res.1
## X-squared = 5.9061, df = 10, p-value = 0.8231
## 
## 
##  Box-Ljung test
## 
## data:  res.1
## X-squared = 13.687, df = 14, p-value = 0.4733
## 
## 
##  Box-Ljung test
## 
## data:  res.2
## X-squared = 4.0006, df = 6, p-value = 0.6766
## 
## 
##  Box-Ljung test
## 
## data:  res.2
## X-squared = 5.7705, df = 10, p-value = 0.8342
## 
## 
##  Box-Ljung test
## 
## data:  res.2
## X-squared = 14.063, df = 14, p-value = 0.445