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!