## 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.

#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
## g_GDP     0.066   -0.027
## g_DGDP    0.051    0.025
## UNRATE   -0.029   -0.029
## 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
##       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")
}