Problem 1

From FREDII data base download the series GDPC1 (quarterly US real GDP). Transform the series in growth rates. Suppose that real GDP growth rates follow an AR(1) with constant such as:

\[ y_t = c + \phi y_{t-1} + \varepsilon_t \]

Where:

\[\varepsilon_t \sim iiid N(0, \sigma^2) \]

(a) Plot the series

#Set key
fredr_set_key("e257a9f359aec337958d438f38801b61")

#Fetch the series 
dta <- fredr(
      series_id = "GDPC1",
      observation_start = as.Date("1947-01-01"),
      observation_end = as.Date("2021-01-01")
        )
names(dta) <- c("date", "series_id", "GDP")

#Calculate growth rates 
dta %<>% mutate(l_GDP = lag(GDP))  
dta %<>% mutate(g_GDP =100*(GDP - l_GDP)/lag(GDP)) 

#Plot the data 
dta %>%  plot_time_series(date, g_GDP, .smooth = FALSE,
                          .line_color="#cb4154",
                          .interactive = FALSE, 
                          .title= "US GDP Growth Rate 1947-2020, quaterly" 
                          )
## Warning: Removed 1 row(s) containing missing values (geom_path).

(b) Estimate AR(1) with constant using OLS

dta %<>% mutate(l_g_GDP = lag(g_GDP))
dta[is.na(dta)] <- 0

X = as.matrix(cbind(1, dta$l_g_GDP))
Y = as.matrix(dta$g_GDP) 

bh = solve(t(X)%*%X)%*%t(X)%*%Y
beta.hat = as.data.frame(cbind(c("Constant", "Phi_1"), bh))
names(beta.hat) = c("Coeff.", "Est")
beta.hat
##     Coeff.               Est
## 1 Constant 0.667992949368724
## 2    Phi_1 0.122406614650672

(c) Estimate and plot the first 10 autocorrelations

For this exercise we are going to compare our user defined function with the built in one. We have to have in mind two things. First, the ACF at lag “h” is the correlation between the series with at “t” and the series with a lead “t+h”. Such as: \[ \rho(h) = Corr(X_{t+h}, X_t) \] The second thing to have in mind is that the built in ACF in R uses zero meaned series. If the series is not zero meaned it substract the mean.

m_g_GDP = as.matrix(dta$g_GDP) 
m_g_GDP = m_g_GDP - mean(m_g_GDP)
n = length(m_g_GDP)                                     # Length of time series
ACF = 0                                              # Starting an empty vector to capture the auto-correlations.
z = sum(m_g_GDP*m_g_GDP) / n                             # Variance of time series.
ACF[1] = z/z                                         # Normalizing first entry to 1.
for(i in 1:24){                                      # Took 24 points to parallel the output of `acf()`
  lag = m_g_GDP[-c(1:i)]                                  # Introducing lags in the stationary ts.
  clipped = m_g_GDP[1:(length(lag))]                      # Compensating by clipping the tail of the ts.
  ACF[i + 1] = (sum(clipped * lag) / n) / z  # Storing each normalized correlation.
}

par(mfrow = c(1,2))
plot(ACF, type="h", main="ACF Manual calculation"); abline(h = 0) # and manual results (right).
abline(h = c(1,-1) * 1.96/sqrt(length(m_g_GDP)), lty = 2, col="red")
acf(m_g_GDP, type = "correlation")

(d) Obtain and plot the coefficients of the Wold representation of the series

mu <- bh[1] / (1 - bh[2])
er <- Y - X%*%bh
s2 <- er%*%t(er)/length(Y)
varu <- s2/(1-bh[2]^2)

th <- bh[2]^seq(0 , 10, by=1)

plot(th, type="b")

\[ \begin{align} y_t&= c + \phi y_{t-1} + \varepsilon_t \\ &= c + \phi(c + \phi y_{t-2} + \varepsilon_{t-1}) + \varepsilon_t \\ &= (1 +\phi)c + \phi^2 y_{t-2} + \varepsilon_t + \phi \varepsilon_{t-1} \\ &= (1 +\phi)c + \phi^2(c + \phi y_{t-3} + \varepsilon_{t-2} ) + \varepsilon_t + \phi \varepsilon_{t-1} \\ &=(1 +\phi + \phi^2)c + \phi^3y_{t-3 } + \varepsilon_t + \phi \varepsilon_{t-1} +\phi^2 \varepsilon_{t-3} \\ &= ... \\ &= (1 + \phi + \phi^2 + ...+ \phi^{t-1})c + \phi^t y_0 + \varepsilon_t + \varepsilon_t + \phi \varepsilon_{t-1} +\phi^2 \varepsilon_{t-2} + .... \phi^{t-1} \varepsilon_{t} \\ \end{align} \]

Problem 2

Now suppose that GDP admits an AR(2) representation

(a) Estimate the parameters with OLS

dta %<>% mutate(l_l_g_GDP = lag(l_g_GDP))
dta[is.na(dta)] <- 0

X = as.matrix(cbind(1, dta$l_g_GDP,dta$l_l_g_GDP))
Y = as.matrix(dta$g_GDP) 

bh = solve(t(X)%*%X)%*%t(X)%*%Y
beta.hat = as.data.frame(cbind(c("Constant", "Phi_1", "Phi_2"), bh))
names(beta.hat) = c("Coeff.", "Est")
beta.hat
##     Coeff.               Est
## 1 Constant 0.593423915293142
## 2    Phi_1  0.10845723813144
## 3    Phi_2  0.11628803141594

(b) Compute the roots of AR polynomial

rootsAR2 <- abs(roots(bh))
rootsAR2
## [1] 0.442675 0.442675
if (rootsAR2[1]==1) {
    print('The process is not stationary')
   } else { print('The process is stationary')
}
## [1] "The process is stationary"
# Check causality of an AR(1) process:


if (rootsAR2[1]<=1) {
    print('The process is not causal')
}else { print('The process is causal')}
## [1] "The process is not causal"

(c) State the conditions under which the process is causal and stationary

If the roots of the autoregressive process are smaller than 1 the process is not causal. In this case, the roots are inside of the unit root circle so they are not causal

(d) Obtain and plot the coefficients of the Wold representation of y_t

We can define the companion form matrix as:

\[ \begin{bmatrix} y_{t} \\ y_{t-1} \end{bmatrix} = \begin{bmatrix} c \\ 0 \end{bmatrix} + \begin{bmatrix} \phi_1 & \phi_2 \\ 1 & 0 \end{bmatrix} \begin{bmatrix} y_{t-1} \\ y_{t-2} \end{bmatrix} + \begin{bmatrix} \varepsilon_t \\ 0 \end{bmatrix} \] Therefore we just iterate in a for loop the companion matrix and we powerit to obtain the Theta coefficiente of the Wold representaiton. As we can see it decasy quickly.

A <- rbind(c(bh[2], bh[3]),c(1,0))

for (j in 1:10) {
    AA<-A^(j-1);
    th[j]<-AA[1,1]
}

plot(th, type="b")

Problem 3 Consider the following MA process

(a) Show that the process is stationary

\[ \begin{align} y_t&= c + \varepsilon_t + 1.2\varepsilon_{t-1} + 2\varepsilon_{t-2} \\ E(y_t)&= E(c + \varepsilon_t + 1.2\varepsilon_{t-1} + 2\varepsilon_{t-2})\\ &=E(c)+E(\varepsilon_{t}) +1.2 E(\varepsilon_{t-1}) + 2E(\varepsilon_{t-2}) = c\\ VAR(y_t)&= E(c + \varepsilon_t + 1.2\varepsilon_{t-1} + 2\varepsilon_{t-2}-c)^2 \\ &= (1+1.2^2 +2^2)\sigma^2\\ \end{align} \]

(b) Is this process invertible?

\[ \begin{align} y_t&= c + \varepsilon_t + 1.2\varepsilon_{t-1} + 2\varepsilon_{t-2} \\ &= (1 + 1.2L +2L^2)\varepsilon_{t}\\ \end{align} \] The roots can be defined as:

\[ \begin{align} 2z^2 + 1.2z + 1 &= 0 \\ z_{1,2} &= \frac{-b \pm \sqrt{b^2 -4ac}}{2a} \\ z_{1,2} &= \frac{-1 \pm \sqrt{1.2^2 -4*2*1}}{2*2} \\ z_{1} &=-0.6+1.280625i \\ z_{2} &= -0.6-1.280625i \end{align} \] We can see that the roots are smaller than 1 in absolute value therefore the process is not invertible.

x1 <- c(2, 1.2, 1) 

# Calling polyroot() function 
polyroot(x1) 
## [1] -0.6+1.280625i -0.6-1.280625i