Problem 1

(1) Create the spread as the difference between the long and the short rate. Create the growth rates of real GDP and the GDP deflator. Plot the series.

First we download, transform and plot the series

#Set key
fredr_set_key("e257a9f359aec337958d438f38801b61")

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

dta_1 <- fredr(
      series_id = c("GDPDEF"),
      observation_start = as.Date("1947-01-01"),
      observation_end = as.Date("2019-01-01")
        )
names(dta_1) <- c("date", "series_id", "DGDP")

dta_2 <- fredr(
      series_id = c("FEDFUNDS"),
      observation_start = as.Date("1947-01-01"),
      observation_end = as.Date("2019-01-01")
        )
names(dta_2) <- c("date", "series_id", "FEDFUNDS")

dta_3 <- fredr(
      series_id = c("UNRATE"),
      observation_start = as.Date("1947-01-01"),
      observation_end = as.Date("2019-01-01")
        )
names(dta_3) <- c("date", "series_id", "UNRATE")

dta_4 <- fredr(
      series_id = c("GS10"),
      observation_start = as.Date("1947-01-01"),
      observation_end = as.Date("2019-01-01")
        )
names(dta_4) <- c("date", "series_id", "GS10")

list(dta,dta_1, dta_2,dta_3, dta_4) %>% reduce(left_join, by = "date") -> data


#Calculate growth rates 
data %<>% mutate(l_GDP = lag(GDP))  
data %<>% mutate(g_GDP =100*(GDP - l_GDP)/lag(GDP)) 
data %<>% mutate(l_DGDP = lag(DGDP))  
data %<>% mutate(g_DGDP =100*(DGDP - l_DGDP)/lag(DGDP)) 
data %<>% mutate(spread = GS10 - FEDFUNDS) 


#Plot the data 
data %>%  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).

data %>%  plot_time_series(date, g_DGDP, .smooth = FALSE,
                          .line_color="#4199cb",
                          .interactive = FALSE, 
                          .title= "US GDP Deflator Growth Rate 1947-2020, quaterly" 
                          )
## Warning: Removed 1 row(s) containing missing values (geom_path).

data %>%  plot_time_series(date, spread, .smooth = FALSE,
                          .line_color="#41cbb8",
                          .interactive = FALSE, 
                          .title= "Spread long and short tame rate , quaterly" 
                          )
## Warning: Removed 30 row(s) containing missing values (geom_path).

#### (2) Using the growth rates of real GDP, the unemployment rate, the inflation rate, the federal funds rate and the spread estimate a VAR(4).

To consider an example, let us assume that there are no deterministic regressors, K= 4 (four variables) and P= 4 (four lags). This can be written as \[ \begin{eqnarray}\nonumber {\bf{y}}_{t}= A_{1} {\bf{y}}_{t-1} + {\bf{u}}_{t} \end{eqnarray} \]

Where the matrices are:

\[ \bf{y}_{t}=\left[ \begin{array}{c} y_{1,t} \\ y_{2,t} \\ y_{3,t} \\ y_{4,t} \\ y_{5,t} \\ \end{array} \right] \] And \[ A_{1}=\left[ \begin{array}{c} \alpha_{11} & \alpha_{12} & \alpha_{12} & \alpha_{13} & \alpha_{14}& \alpha_{15} \\ \alpha_{11} & \alpha_{22} & \alpha_{22} & \alpha_{23} & \alpha_{24}& \alpha_{25} \\ \alpha_{11} & \alpha_{32} & \alpha_{22} & \alpha_{33} & \alpha_{34}& \alpha_{35} \\ \alpha_{11} & \alpha_{42} & \alpha_{42} & \alpha_{43} & \alpha_{44}& \alpha_{45} \\ \alpha_{11} & \alpha_{52} & \alpha_{52} & \alpha_{53} & \alpha_{54}& \alpha_{55} \\ \end{array} \right] \] And the error \[ \bf{u}_{t}=\left[ \begin{array}{c} u_{1,t} \\ u_{2,t} \\ u_{3,t} \\ u_{4,t} \\ u_{5,t} \\ \end{array} \right] \]

data %>% select(g_GDP, g_DGDP, UNRATE, spread, FEDFUNDS)-> data_var
data_var %<>% slice(31:nrow(data_var))

p <- 4                          #Define the number of lags
T <- as.numeric(nrow(data_var)) #Define the sample size

p_1 <- p +1
Y = data.matrix(slice(data_var, p_1:T)) #Create Y matrix

data_var_X <- as.data.frame(1)  #Creat the row vector with the 
names(data_var_X) <- "const"    #Constant

for(i in 1:p) {
    T <- T-1
    data_var_X <- cbind(data_var_X, slice(data_var, p:T))
    p <- p-1
}

X = data.matrix(data_var_X)      #Create X matrix

Once we have defined the matrices, we can estimate the VAR using the OLS estimator.

p <- 4                          #Define the number of lags
T <- nrow(data_var)
N <- as.numeric(ncol(data_var))
phi_hat = round(solve(t(X)%*%X)%*%t(X)%*%Y, digits = 3)
Y_hat = X%*%phi_hat
err = Y - Y_hat
omega = t(err)%*%err/(T-N*p-1)
phi_hat
##           g_GDP g_DGDP UNRATE spread FEDFUNDS
## const    -0.085  0.038  0.486 -0.175   -0.007
## g_GDP     0.218  0.018 -0.174 -0.134    0.342
## g_DGDP    0.266  0.464 -0.086 -0.318    0.484
## UNRATE    0.018 -0.109  0.986  0.647   -0.527
## spread    0.059  0.001  0.008  0.468    0.467
## FEDFUNDS -0.115  0.033  0.025 -0.105    1.056
## g_GDP     0.293 -0.001 -0.133 -0.091    0.105
## g_DGDP    0.253  0.210 -0.046 -0.031    0.049
## UNRATE    0.506  0.250  0.006 -0.763    0.629
## spread   -0.062 -0.011 -0.019  0.292   -0.380
## FEDFUNDS -0.046 -0.012 -0.024  0.147   -0.217
## g_GDP     0.096  0.021 -0.004  0.019   -0.023
## g_DGDP   -0.429  0.079  0.068 -0.307    0.491
## UNRATE   -0.244 -0.121 -0.025  0.520   -0.555
## spread    0.055 -0.057  0.034  0.069   -0.037
## FEDFUNDS  0.193 -0.054  0.045 -0.076    0.163
## g_GDP     0.090  0.063 -0.012  0.135   -0.098
## g_DGDP   -0.074  0.209  0.056  0.186   -0.341
## UNRATE   -0.255 -0.024 -0.019 -0.219    0.334
## spread    0.066  0.051 -0.029 -0.189    0.118
## FEDFUNDS -0.027  0.025 -0.029  0.011   -0.054
omega
##                g_GDP       g_DGDP       UNRATE      spread    FEDFUNDS
## g_GDP     0.58432783 -0.011876787 -0.097081243  0.01036834  0.03740627
## g_DGDP   -0.01187679  0.065288026 -0.001241645 -0.03558327  0.05346912
## UNRATE   -0.09708124 -0.001241645  0.078711156  0.04008611 -0.06291717
## spread    0.01036834 -0.035583270  0.040086106  0.78053680 -0.73996980
## FEDFUNDS  0.03740627  0.053469121 -0.062917167 -0.73996980  0.92061923

(3) Compute and plot the impulse response functions of the Wold shocks (the Wold representation

The companion-form It is possible to express a VAR(p) model in a VAR(1) form, which turns out to be very useful in many practical situations (and for some technical derivations). This reformulation is done by expressing the VAR(p) model in the companion-form, which is accomplished as follows. Consider the model:

\[ \begin{eqnarray}\nonumber Z_{t}=\Gamma_{0}+\Gamma_{1}Z_{t-1}+ \Upsilon_{t} \end{eqnarray} \] Where we have defined:

$$ \[\begin{eqnarray}\nonumber Z_{t}=\left[ \begin{array}{c} {\bf{y}}_{t} \\ {\bf{y}}_{t-1} \\ \vdots \\ {\bf{y}}_{t-p+1} \end{array} \right] , \hspace{1cm} \Gamma_0=\left[ \begin{array}{c} \mu \\ 0 \\ \vdots \\ 0 \end{array} \right] , \hspace{1cm} \Upsilon_{t} =\left[ \begin{array}{c} {\bf{u}}_{t} \\ 0 \\ \vdots \\ 0 \end{array} \right] \end{eqnarray}\]

\[ Which can be translated as matrix notation as: \] \[\begin{eqnarray}\nonumber \left[ \begin{array}{c} {\bf{y}}_{t} \\ {\bf{y}}_{t-1} \\ {\bf{y}}_{t-2} \\ \vdots \\ {\bf{y}}_{t-p+1} \end{array} \right] =\left[ \begin{array}{c} \mu \\ 0 \\ 0 \\ \vdots \\ 0 \end{array} \right] +\left[ \begin{array}{ccccccc} A_{1} & A_{2} & \cdots & A_{p-1} & A_{p} \\ I & 0 & \cdots & 0 & 0 \\ 0 & I & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \cdots & I & 0 \end{array} \right] \left[ \begin{array}{c} {\bf{y}}_{t-1} \\ {\bf{y}}_{t-2} \\ {\bf{y}}_{t-3} \\ \vdots \\ {\bf{y}}_{t-p} \end{array} \right] + \left[ \begin{array}{c} {\bf{u}}_{t} \\ 0 \\ 0 \\ \vdots \\ 0 \end{array} \right] \end{eqnarray}\] $$

c = 1

phi_hat_aux = t(phi_hat[2:nrow(phi_hat),])
m_Diagonal = diag(N*p-N)
m_Zeros  = matrix(0,N*p-N, N)

Capital_A <- rbind(phi_hat_aux, cbind(m_Diagonal, m_Zeros))
head(Capital_A)
##           g_GDP g_DGDP UNRATE spread FEDFUNDS  g_GDP g_DGDP UNRATE spread
## g_GDP     0.218  0.266  0.018  0.059   -0.115  0.293  0.253  0.506 -0.062
## g_DGDP    0.018  0.464 -0.109  0.001    0.033 -0.001  0.210  0.250 -0.011
## UNRATE   -0.174 -0.086  0.986  0.008    0.025 -0.133 -0.046  0.006 -0.019
## spread   -0.134 -0.318  0.647  0.468   -0.105 -0.091 -0.031 -0.763  0.292
## FEDFUNDS  0.342  0.484 -0.527  0.467    1.056  0.105  0.049  0.629 -0.380
##           1.000  0.000  0.000  0.000    0.000  0.000  0.000  0.000  0.000
##          FEDFUNDS  g_GDP g_DGDP UNRATE spread FEDFUNDS  g_GDP g_DGDP UNRATE
## g_GDP      -0.046  0.096 -0.429 -0.244  0.055    0.193  0.090 -0.074 -0.255
## g_DGDP     -0.012  0.021  0.079 -0.121 -0.057   -0.054  0.063  0.209 -0.024
## UNRATE     -0.024 -0.004  0.068 -0.025  0.034    0.045 -0.012  0.056 -0.019
## spread      0.147  0.019 -0.307  0.520  0.069   -0.076  0.135  0.186 -0.219
## FEDFUNDS   -0.217 -0.023  0.491 -0.555 -0.037    0.163 -0.098 -0.341  0.334
##             0.000  0.000  0.000  0.000  0.000    0.000  0.000  0.000  0.000
##          spread FEDFUNDS
## g_GDP     0.066   -0.027
## g_DGDP    0.051    0.025
## UNRATE   -0.029   -0.029
## spread   -0.189    0.011
## FEDFUNDS  0.118   -0.054
##           0.000    0.000
tail(Capital_A)
##  g_GDP g_DGDP UNRATE spread FEDFUNDS g_GDP g_DGDP UNRATE spread FEDFUNDS g_GDP
##      0      0      0      0        0     0      0      0      0        1     0
##      0      0      0      0        0     0      0      0      0        0     1
##      0      0      0      0        0     0      0      0      0        0     0
##      0      0      0      0        0     0      0      0      0        0     0
##      0      0      0      0        0     0      0      0      0        0     0
##      0      0      0      0        0     0      0      0      0        0     0
##  g_DGDP UNRATE spread FEDFUNDS g_GDP g_DGDP UNRATE spread FEDFUNDS
##       0      0      0        0     0      0      0      0        0
##       0      0      0        0     0      0      0      0        0
##       1      0      0        0     0      0      0      0        0
##       0      1      0        0     0      0      0      0        0
##       0      0      1        0     0      0      0      0        0
##       0      0      0        1     0      0      0      0        0

Once we have the the Companion Form, we check the stability and sstationaryti.

abs_eigen_Capital_A <- abs(eigen(Capital_A)$values)

if (max(abs_eigen_Capital_A)>1){
print('There is at least 1 eigenvalue of the companion form matrix larger than 1, hence the VAR is not stable.')
} else {
print('The largest eigenvalue is smaller than 1, so the VAR is stable and stationary.')
}
## [1] "The largest eigenvalue is smaller than 1, so the VAR is stable and stationary."

Now we do the Wold representation impulse responses

horizon=20; # horizons j=1,...,20
someData <- rep(0, N*N*horizon)
C <- array(someData, c(N, N, horizon)) # tensor: 3 dimensional object that stores n*n matrices for j=0,1,...,20

for (j in 1:horizon){
    Capital_C=Capital_A^(j-1)
    C[,,j]=Capital_C[1:N,1:N] 
}

Oganize hte matrices in a horizon x N*N matrix:

Intermediate_C <- aperm(C, c(3,2,1))
C_Wold <-Reshape(Intermediate_C, horizon, N*N)

Now we can plot the IRF

for(i in 1:ncol(C_Wold)){
  plot(seq(1:horizon), C_Wold[,i], type="l")
}