[R] nlminb error

[Rodrigo Andres Cristancho Castellanos]

Hi, I´m having troubel using nlminb this is the warning that shows up.

Warning: Cholesky factorization 'dpotrf' exited with status: 1
Variance of the prediction error can not be computed.
Warning: Cholesky factorization 'dpotrf' exited with status: 1
Determinant of the variance of the prediction error can not be computed.

I´m trying to optimize a likelihood function. The code is a little messy,
you can find it at the end of the email.

I´ll appreciate any kind of help.

Regards.

Rodrigo Cristancho.



library(foreach)
library(FKF)
library(ggplot2)
library(gridExtra)

# Model setup
#3 Factor Independent Model

delta1<- c(0.006, 0.003, 0.007)
kappa <- c(0.001, 0.002, 0.004)
sigma <- c(0.002, 0.005, 0.003)
rc    <- c(0.00002)
r1    <- c(0.00002)
r2    <- c(0.5)

#All other parameters
T         <- 110  # years
delta_t   <- 1 # monthly zeros observations
dt        <- 1 # monthly r simulations
n         <- T/dt # number of r simulations
r_0       <- delta1
measurement_error <- 0.001 # for zero-coupon rates
m         <- length(delta1)   # dimension of state variables
# maturity  <- c(1/12 ,1/4, 1/2, 10) # zeros for 1 factor simulation
# maturity  <- c(1/12 ,1/4, 1/2, 1, 2, 3, 5, 7, 10) # zeros for 2 factor
simulation

maturity  <- 61-xc #
d         <- length(maturity)   # dimension of observations
N         <- T/delta_t # number of observations
premubar1=list(premubar)

# PARAMETER ESTIMATION
# ----------------------------------------------------------
# Random initialization of parameters
init_params_if <- function()
{
  delta1_init <<- runif(m, min=0.0, max=0.05)
  kappa_init <<- runif(m, min=0.0, max=0.05)
  sigma_init <<- runif(m, min=0.0, max=0.05)
  rc_init <<- runif(1, min=0.0, max=0.001)
  r1_init <<- runif(1, min=0.0, max=0.001)
  r2_init <<- runif(1, min=0.0, max=0.001)
}
# optimization parameter bounds

upper_bound <- c(rep(c(1.0, 1.0, 1.0, 1.0),each=m), rep(0.1,d))
lower_bound <- c(rep(c(0.0001, 0.0001, 0.0001,-1.0 ),each=m), rep(0.0001,d))
actual_param <- c(kappa=kappa, delta1=delta1, sigma=sigma, rc=rep(rc,1),
r1=rep(r1,1), r2=rep(r2,1))

#
#------------------------------------------------------------------------------------------------------------
#


#Kalman Filter for 3 Factor Indipendent Model
Ind_Fact_KF <- function(delta1, kappa, sigma, rc, r1, r2, observations)
{
  # initial state variable (a0: m x 1)
  #r_init <- as.vector(delta1)                    # unconditional mean of
state variable

  r_init <- as.vector(c(rep(0,m)))

  # variance of state variable (P0: m x m)
  #P_init <- (sigma^2/(2*(delta1^3)))*diag(1,m,m)   # unconditional
variance of state variable

  P_init <- (sigma^2/(2*(delta1)))*diag(1,m,m)

  # intercept of state transition equation (dt: m x 1)
  C <- matrix(0,nrow = m,ncol = 1)

  # factor of transition equation (Tt: m x m x 1)
  F_ <- array(exp(-kappa*delta_t)*diag(m),dim=c(m,m,1))

  # factor of measurement equation (Zt: d x m x 1)
  B <-
array(1/matrix(rep(delta1,d),d,byrow=TRUE)*(1-exp(-matrix(rep(delta1,d),d,byrow=TRUE)
* matrix(rep(maturity,m),d))),dim=c(d,m,1))

  # intercept of measurement equation (ct: d x 1)

  A <-
(1/2)*t(sigma^2/delta1^3)%*%t(((1/2)*(1-exp(-2*(matrix(rep(delta1,d),d,byrow=TRUE)*matrix(rep(maturity,m),d)))))

 -2*(1-exp(-2*(matrix(rep(delta1,d),d,byrow=TRUE)*matrix(rep(maturity,m),d))))
          -(matrix(rep(delta1,d),d,byrow=TRUE)*matrix(rep(maturity,m),d)))

  A <- matrix(-t(A)/maturity,nrow=length(maturity))

  B <- array(B[,,1]/matrix(rep(maturity,m),d),dim=c(d,m,1))

  # variance of innovations of transition (HHt: m x m x 1)
  Q <-
array(sigma^2/(2*kappa)*(1-exp(-2*kappa*delta_t))*diag(m),dim=c(m,m,1))

  # variance of measurement error (GGt: d x d x 1)

  R <- array(diag(d)*(rc+r1*exp(r2)),dim=c(d,d,1))
  #R <- array(diag(d)*(rc),dim=c(d,d,1))

  ##Funcion de filtro de Kalman rapido con base al desarrollo del modelo
arriba
  filtered_process <- fkf(a0=r_init, P0=P_init, dt=C, ct=A, Tt=F_, Zt=B,
HHt=Q, GGt=R, yt=observations)
  return(filtered_process)
}

#aaaa=fkf(a0=r_i(nit, P0=P_init, dt=C, ct=A, Tt=F_, Zt=B, HHt=Q, GGt=R,
yt=t(premubar))
aaaa=Ind_Fact_KF(delta1, kappa, sigma, rc, r1, r2,observations=t(premubar))
aaaa$logLik
aaaa$Ft

# Retrieve short rates using Kalman Filter
retrieve_short_rates_if <- function(rates, optim_controls,
lower_bound=NULL, upper_bound=NULL)
{
  observations <- rates
  init_params_if()
  initial_param <<- c(delta1=delta1_init,kappa=kappa_init,
sigma=sigma_init, rc=rc_init, r1=r1_init, r2=r2_init)

  if_KF_loglik <- function(x)
  {
    delta1 <- x[1:m]; kappa <- x[(m+1):(2*m)]; sigma <- x[(2*m+1):(3*m)];
rc <- x[(3*m+1):(3*m+1)]; r1 <- x[(3*m+2):(3*m+2)];r2 <-
x[(3*m+3):length(x)]
    return(-Ind_Fact_KF(delta1,kappa,sigma,rc,r1,r2,observations)$logLik)
  }

  # optimization of log likelihood function
  fitted_model <- nlminb(initial_param, if_KF_loglik,
control=optim_controls, lower=lower_bound, upper=upper_bound)
  return(fitted_model)
}

rrif=retrieve_short_rates_if(rates=t(premubar),optim_controls=optim_controls,upper=upper_bound,lower=lower_bound)

	[[alternative HTML version deleted]]