I am trying to estimate coefficients of a state-space model described in [Diebold et.al (2006)][1] with data and scripts [here][2]:
$$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][3] with a mapping function [here][4]. The guide is quite in-depth and extensive, however the implementation in MATLAM 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][5]:
$$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 <- vecm
      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!


  [1]: https://www.sas.upenn.edu/~fdiebold/papers/paper55/DRAfinal.pdf
  [2]: http://www.estima.com/procs_perl/dra_joe_2006.zip
  [3]: https://de.mathworks.com/help/econ/using-the-kalman-filter-to-estimate-and-forecast-the-diebold-li-model.html
  [4]: https://www.mathworks.com/matlabcentral/mlc-downloads/downloads/submissions/47479/versions/2/previews/Example_DieboldLi.m/index.html
  [5]: https://stats.stackexchange.com/questions/21554/dlm-package-state-equation-maximum-likelihood-with-constant-terms