# [R] nlminb error

Rodrigo Andres Cristancho Castellanos cristanchocc at gmail.com
Sun Apr 9 02:21:00 CEST 2017

```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]]

```