# MLE estimation of mean-adjusted state-space model

I am trying to estimate coefficients of a state-space model described in Diebold et.al (2006) with data and scripts here: $$y_t = Zf_t + \epsilon_t$$ $$f_t-\mu_t = T(f_{t-1}-\mu)+\eta_t$$ The main issue is that I have not been able to specify this model in the dlm package for R, as the estimation of the means is more complicated than expected. I have found a guide for MATLAB with a mapping function here. The guide is quite in-depth and extensive, however the implementation in MATLAB and dlm package differ. The two-step VAR(1) approach is the same:

rm(list=ls())
setwd(dirname(rstudioapi::getSourceEditorContext()$path)) require(dlm) df <- read.delim('dra data.txt') yield <- as.matrix(df[,c(2:(ncol(df)-3))]) maturities <- c(3,6,9,12,15,18,21,24,30,36,48,60,72,84,96,108,120) lambda0 <- 0.0609 X <- matrix(c(rep(1,length(maturities)), (1-exp(-lambda0*maturities))/(lambda0*maturities), ((1-exp(-lambda0*maturities))/(lambda0*maturities) - exp(-lambda0*maturities))), ncol = 3) beta <- matrix(rep(0,3*nrow(yield)), ncol = 3) eps <- matrix(rep(0,length(maturities)*nrow(yield)), ncol = length(maturities)) i <-1 for (i in 1:nrow(yield)) { y <- as.data.frame(yield[i,]) data <- cbind(y,X) names(data)<- c("y", "c1","c2","c3") reg <- lm(formula = y~c1+c2+c3-1, data = data) beta[i,]<- reg$$coefficients eps[i,] <- reg$$residuals } colnames(beta)<- c("b0","b1","b2") library(vars) VAR(beta,p = 1, type = "const")->varbeta varbeta$varresult
varlagbeta <- matrix(c(varbeta[["varresult"]][["b0"]][["coefficients"]][1:3],
varbeta[["varresult"]][["b1"]][["coefficients"]][1:3],
varbeta[["varresult"]][["b2"]][["coefficients"]][1:3]),byrow = T, ncol = 3)
rownames(varlagbeta)<-c("beta_0","beta_1","beta_2")
colnames(varlagbeta)<-c("beta_0_l1","beta_1_l1","beta_2_l1")

cat('VAR(1) matrix of estimateed coefficients: \n')
print(varlagbeta)
vecG <- as.vector(t(varlagbeta))
cat("VAR(1) covariance matrix of residuals: \n")
print(summary(varbeta)$covres)  The one-step state-space model is not so successful. I've rewritten the model as in the tutorial and this post: $$y_t-Z\mu = Zx_t+e_t$$ $$x_t = Tx_{t-1}+\eta_t$$ Where $$x_t = f_{t-1}-\mu$$ and tried to implement it in dlm: vecG <- as.vector(t(varlagbeta)) matW <- summary(varbeta)$covres
matW <- chol(matW)
vecW <- c(log(matW[1,1]),matW[1,2],matW[1,3],log(matW[2,2]),matW[2,3],log(matW[3,3]))
vecV <- diag(cov(eps))
vecC <- cov(beta)
vecm <- apply(beta,2,mean)

params0 <- c(vecG, vecW, log(vecV), vecm, log(lambda0))
param<- params0
var_rest <- function(x){
return(exp(x))
}

stl <- ncol(yield)
poz <- nrow(yield)

estpar <- function(param){
lam <- var_rest(param[length(param)])

F.mat <-matrix(rep(0,3*stl),nr=stl)
F.mat[,1:3] <- c(rep(1,length(maturities)), (1-exp(-lam*maturities))/(lam*maturities),
((1-exp(-lam*maturities))/(lam*maturities) - exp(-lam*maturities)))

V <- diag(var_rest(param[16:32]))

G.mat <- matrix(rep(0,3*3),nr=3)

G.mat[1:3, 1:3] <- matrix(param[1:9],nrow =3, ncol = 3, byrow = TRUE)

W.mat <-matrix(rep(0,3*3),nrow=3)
#param[c(28,31,33)] <- var_rest(param[c(28,31,33)])
param[c(10,13,15)] <- var_rest(param[c(10,13,15)])
W.mat[1,1] <-param[c(10)]
W.mat[2,1:2] <-param[c(11,13)]
W.mat[3,1:3] <-param[c(12,14,15)]
W <- W.mat%*%t(W.mat)

m0.mat <- rep(0,3)
C0.mat <- diag(1e6,3)

return( dlm(m0=m0.mat, C0=C0.mat, FF=F.mat, GG=G.mat, W=W,V=V))
}

dlm_optim <- function(y, parm, build, method = "BFGS", ..., debug = FALSE){
logLik <- function(parm,y, ...) {
mod <- build(parm, ...)
y2 <- sweep(y,2,X%*%parm[33:35])
return(dlmLL(y = y2, mod = mod, debug = debug))
}
out <- optim(parm, logLik, y=y, method = method, ...)
return(out)
}

pok1a <- dlm_optim(as.matrix(yield),params0,build = estpar, hessian=T,control=list(maxit=10000))#, control=list(maxit=5))


MMy estimates are completely off. There are 3 main issues:

1. How to specify the estpar function so that the means will be estimated as well?
2. What to do with the $$m_0$$ and $$C_0$$ parameters? The Matlab example does not specify anything like it, but the dlm function requires it as input.
3. If I use BFGS optimization method (which the authors used), I have singularity issues with the covariance matrices. If I use L-BFGS-B, the results are off.

I know this is a long and complicated question, but any help would be appreaciated!

• I've found that there is much bigger issue with the extended model they use - the authors incorporate macro-variables as observed state variables! So these variables (capacity utilization, inflation and FFR, using your notation) are in the observation vector $y_t$, as well as in the state vector $\alpha_t$. This I could get around by extending matrix $Z$, however, then the matrix $H$ would not be invertible (as the errors of ! I think I will create a new question with this and share the link in edit. Commented Dec 29, 2020 at 18:08