I am trying to estimate the following state space model:
\begin{equation} y_t = y^{gap}_{t} + y^*_t \end{equation} \begin{equation} y^{gap}_{t} = \alpha_{1}y^{gap}_{t-1}+\alpha_{3}y^{gap}_{t-2} +\alpha_{2}/2(r_{t-1}-r^*_{t-1}) + \alpha_{2}/2(r_{t-2}-r^*_{t-2}) + \epsilon^{ygap}_t \end{equation} \begin{equation} y^*_t = y^*_{t-1}+ \mu_{t-1} + \epsilon^{y^*}_t \end{equation} \begin{equation} \mu_t = \mu_{t-1} + \epsilon^{\mu}_t \end{equation} \begin{equation} u_t = u^{gap}_{t} + u^*_t \end{equation} \begin{equation} u^{gap}_{t} = \gamma_10.4y^{gap}_{t}+\gamma_10.3y^{gap}_{t-1}+\gamma_10.2y^{gap}_{t-2}+\gamma_10.1y^{gap}_{t-3} +\epsilon^{u^{gap}}_t \end{equation} \begin{equation} u^*_t = u^*_{t-1} + \epsilon^{u^*}_t \end{equation} \begin{equation} \pi_{t} = \beta_{1}/3\pi_{t-1} + \beta_{1}/3\pi_{t-2} + \beta_{1}/3\pi_{t-3} + \beta_{2}u^{gap}_{t-1} + (1-\beta_{1})\pi^{e}_{t} +\epsilon^{\pi}_t \end{equation} \begin{equation} r^*_t = 4\mu_t+z_{t} \end{equation} \begin{equation} z_t = z_{t-1} +\epsilon^{z}_t \end{equation}
However, the results I am getting are completely wrong. For example, my output gap is trending like a non-stationary variable.
I have had a crack a casting this model into a form suitable for R's DLM package:
\begin{equation} Y_t = F\theta_t + V_t \end{equation} \begin{equation} \theta_t = G\theta_{t-1} + W_t \end{equation}
Where
\begin{equation} Y = \left[\begin{array}{c} y_{t} \\ u_{t} \\ r_{t} \\ \pi_{t} \end{array}\right] \end{equation}
\begin{equation} F = \left[\begin{array}{ccccccccccccccccc} 1 & 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 & 1 & 1 & 0 & 0 & 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 & 1 & 0 & 0 & 0\\ \end{array}\right] \end{equation}
\begin{equation} \theta_{t} = \left[\begin{array}{c} y^* \\y^{gap}_{t} \\y^{gap}_{t-1} \\y^{gap}_{t-2} \\ \mu_{t} \\z_{t}\\ r^*_{t} \\ r^*_{t-1}\\ r_{t} \\ r_{t-1} \\ u^*_{t}\\ u^{gap}_{t}\\ u^{gap}_{t-1} \\ \pi_{t} \\\pi_{t-1} \\ \pi_{t-3} \\ 1-\beta_1 \end{array} \right] \end{equation}
\begin{equation} G =\left[\begin{array}{ccccccccccccccccc} 1& 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & \alpha_1 & \alpha_2 & 0 & 0 & 0 & \alpha_3/2 & \alpha_3/2 & \alpha_3/2 & \alpha_3/2 & 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 & 1 & 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 & 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 & 4 & 1 & 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 & 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 & 1 & 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 & \gamma_1(0.4\alpha_1+0.3) & \gamma_1(\alpha_20.4+0.2) & \gamma_10.1 & 0 & 0 & \gamma_10.4\alpha_3/2 & \gamma_10.4\alpha_3/2 & \gamma_10.4\alpha_3/2 & \gamma_10.4\alpha_3/2 & 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 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \beta_3 & \beta_1/3 & \beta_1/3 & \beta_1/3 & \pi^{e}_{t} \\ 0 & 0 & 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 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\ \end{array}\right] \end{equation} \begin{equation} R =\left[\begin{array}{ccccccccccccccccc} \sigma^{y^*} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & \sigma^{y^{gap}} & 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 & \sigma^{\mu} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & \sigma^{z} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 4\sigma^{\mu} & \sigma^{z} & 4\sigma^{\mu}+\sigma^{z} & 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 & \sigma^{r_t} & 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 & \sigma^{u^{*}_t} & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & \gamma_1\sigma^{y^{gap}_t} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \sigma^{u^{gap}_t}+\gamma_1\sigma^{y^{gap}_t} & 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 & \sigma^{\pi} & 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\\ \end{array}\right] \end{equation} \begin{equation} W=RR^T \end{equation}
As the above shows, I am using a little trick by treating r as a random walk - as advised by the author of the DLM package in response to a question: https://r.789695.n4.nabble.com/Setting-up-a-State-Space-Model-in-dlm-td3580664.html
My code is below, I have used the estimates from this paper for starting values of the coefficients: https://www.rba.gov.au/publications/bulletin/2017/sep/pdf/bu-0917-2-the-neutral-interest-rate.pdf
I apologise as this is not 100% reproducable as I am not sure how to post the data for the model - I am happy to share! Code:
NRDLM <- dlm(
FF = matrix(0,4,17),
V = diag(0.00001, 4),
GG =diag(0,17),
JGG = diag(1,17),
W = diag(0,17),
m0 = rep(0,17),
C0 = diag(1000000,17),
X = NRdata[,c("Inflation.e")]
)
# Matrix to parametrise VCV matrix W
R <- diag(0,17)
# Set all elements of JGG to zero (will change below)
NRDLM$JGG <- diag(0,17)
# build DLM
buildNRDLM <- function(p){
FF(NRDLM)[1,1] <- 1
FF(NRDLM)[1,2] <- 1
FF(NRDLM)[2,11] <- 1
FF(NRDLM)[2,12] <- 1
FF(NRDLM)[3,9] <- 1
FF(NRDLM)[4,14] <- 1
GG(NRDLM)[1,1] <- 1
GG(NRDLM)[1,5] <- 1
GG(NRDLM)[2,2] <- p[1]
GG(NRDLM)[2,3] <- p[2]
GG(NRDLM)[2,7] <- p[3]/2
GG(NRDLM)[2,8] <- p[3]/2
GG(NRDLM)[2,9] <- -p[3]/2
GG(NRDLM)[2,10] <- -p[3]/2
GG(NRDLM)[3,2] <- 1
GG(NRDLM)[4,3] <- 1
GG(NRDLM)[5,5] <- 1
GG(NRDLM)[6,6] <- 1
GG(NRDLM)[7,5] <- 4
GG(NRDLM)[7,6] <- 1
GG(NRDLM)[8,7] <- 1
GG(NRDLM)[9,9] <- 1
GG(NRDLM)[10,9] <- 1
GG(NRDLM)[11,11] <- 1
GG(NRDLM)[12,2] <- p[4]*(0.4*p[1]+0.3)
GG(NRDLM)[12,3] <- p[4]*(0.4*p[2]+0.2)
GG(NRDLM)[12,4] <- p[4]*0.1
GG(NRDLM)[12,7] <- p[4]*0.4*p[3]/2
GG(NRDLM)[12,8] <- p[4]*0.4*p[3]/2
GG(NRDLM)[12,9] <- p[4]*0.4*-p[3]/2
GG(NRDLM)[12,10] <- p[4]*0.4*-p[3]/2
GG(NRDLM)[13,12] <- 1
GG(NRDLM)[14,13] <- p[5]
GG(NRDLM)[14,14] <- p[6]/3
GG(NRDLM)[14,15] <- p[6]/3
GG(NRDLM)[14,16] <- p[6]/3
GG(NRDLM)[15,14] <- 1
GG(NRDLM)[16,15] <- 1
GG(NRDLM)[14,17] <- 1
JGG(NRDLM)[14,17] <- 1
# Variance covariance - RR'
R[1,1] <- p[7]
R[2,2] <- p[8]
R[5,5] <- p[9]
R[6,6] <- p[10]
R[6,7] <- p[10]
R[7,5] <- 4*p[9]
R[9,9] <- p[11]
R[11,11] <- p[12]
R[12,12] <- p[13]
R[12,2] <- p[4]*p[8]
R[14,14] <- p[14]
W(NRDLM) <- R%*%t(R)
m0(NRDLM) <- c(NRdata$log.output[1],0,0,0,mean(diff(NRdata$log.output[1:4])),0,NRdata$real.r[2],NRdata$real.r[1],NRdata$real.r[2],NRdata$real.r[1],NRdata$unr[1],0,0,NRdata$Inflation[3],NRdata$Inflation[2],NRdata$Inflation[1],0)
return(NRDLM)
}
theta <- c(1.53,-0.54, -0.05, 0.62, -0.32, 0.39, 0.38, 0.54, 0.05, 0.22, 0 , 0.15, 0.07, 0.79 ) # estimates from paper
# Estimate model
NRDLM.est <- dlmMLE(y = cbind(NRdata$log.output,NRdata$unr,NRdata$real.r,NRdata$Inflation), parm = theta, build = buildNRDLM, lower =c(rep(-Inf,6),rep(exp(-8),7)), upper= c(rep(Inf,6),rep(exp(12),11)),
control = list(trace = 1, REPORT = 5, maxit = 1000), hessian = TRUE, debug = F)
#--------------------------------------------------------------------------------------------------------------------------
# filtered and smoothed estimates
#--------------------------------------------------------------------------------------------------------------------------
NRDLMbuilt <- buildNRDLM(NRDLM.est$par)
filtered <- dlmFilter(y =cbind(NRdata$log.output,NRdata$unr,NRdata$real.r,NRdata$Inflation), mod = NRDLMbuilt)
smoothed <- dlmSmooth(y = cbind(NRdata$log.output,NRdata$unr,NRdata$real.r,NRdata$Inflation), mod = NRDLMbuilt)
Anything obviously wrong here?