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]]
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
