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]]
I suspect, as you hinted, there's little to no hope that anyone will be willing or able to navigate your code. **Usually** (whatever that means!) these sorts of problems can be traced back to overparameterization -- too few data, which could also mean a lot of "correlated" data, chasing too many parameters (leading to ridges in the surface or problems near the boundaries). This is the bane of numerical optimizers. Try fitting simpler models with fewer parameters if possible, at least to see if convergence can be achieved. Better starting values, other optimizers, and/or use of analytical gradients and hessians can also sometimes help. Sorry I can't be more specific. Perhaps someone more knowledgeable than I can be. Cheers, Bert Bert Gunter "The trouble with having an open mind is that people keep coming along and sticking things into it." -- Opus (aka Berkeley Breathed in his "Bloom County" comic strip ) On Sat, Apr 8, 2017 at 5:21 PM, Rodrigo Andres Cristancho Castellanos <cristanchocc at gmail.com> wrote:> 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 at 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.
Bert has already suggested that your code is too spaghetti to let us follow it. I tried and got a msg that PKF isn't available for R 3.3.3. Do you actually have a test that the function can be evaluated? That's often the main issue. And it does help to try a few parameter sets to see if you have a reasonable function. JN On 2017-04-08 08:21 PM, Rodrigo Andres Cristancho Castellanos wrote:> 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 at 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. >
Whether you have converged to a global optimum or not depends on the nature of the likelihood surface. Have you tried different starting values? As for the warning, I leave that to Prof. Nash, as he *is* (way) more knowledgeable. However, I suspect the answer is yes, it is concerning. Are you at a boundary? Cheers, Bert Bert Gunter "The trouble with having an open mind is that people keep coming along and sticking things into it." -- Opus (aka Berkeley Breathed in his "Bloom County" comic strip ) On Sun, Apr 9, 2017 at 5:55 PM, Rodrigo Andres Cristancho Castellanos <cristanchocc at gmail.com> wrote:> Hi Bert and JC, thanks for the quick response. > > Thinking about what Bert mentions as posible "issues" I think the three of > them are feasible explanations of the error. > > Anyway, I believe I have already found a work around, I just replace the > function "nlminb" with the function "optim". Even though, the warning > persits the result seems to converge without other issues, given that it > indicates that the convergence is succesfully completed. > > The only question that comes to me with this lucky strike, is, do I need > keep worring about the warning ? > > Thank you both for your help. > > Regards, > > Rodrigo Cristancho > > > > 2017-04-09 12:18 GMT-05:00 Bert Gunter <bgunter.4567 at gmail.com>: >> >> I suspect, as you hinted, there's little to no hope that anyone will >> be willing or able to navigate your code. **Usually** (whatever that >> means!) these sorts of problems can be traced back to >> overparameterization -- too few data, which could also mean a lot of >> "correlated" data, chasing too many parameters (leading to ridges in >> the surface or problems near the boundaries). This is the bane of >> numerical optimizers. >> >> Try fitting simpler models with fewer parameters if possible, at least >> to see if convergence can be achieved. Better starting values, other >> optimizers, and/or use of analytical gradients and hessians can also >> sometimes help. >> >> Sorry I can't be more specific. Perhaps someone more knowledgeable >> than I can be. >> >> Cheers, >> Bert >> >> >> Bert Gunter >> >> "The trouble with having an open mind is that people keep coming along >> and sticking things into it." >> -- Opus (aka Berkeley Breathed in his "Bloom County" comic strip ) >> >> >> On Sat, Apr 8, 2017 at 5:21 PM, Rodrigo Andres Cristancho Castellanos >> <cristanchocc at gmail.com> wrote: >> > 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 at 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. > >