Hello,
I have been judiciously using RCPP for six months, and have had my numerous
questions answered so far from previous threads. I have now run into an
issue I cannot get a handle on:
I have converted a fairly large (recursive) Kalman filter function from R
to C++ using RCPP. This function is then called within the R optim()
function, subject to constraints, which estimates a set of ~30 parameters.
My issue occurs within the C++ Kalman Filter function, where, I use the
optimize() function to calculate the yield to maturity rate of various
bonds. I do this by calling back to the R environment to access the
optimize() function. Below is the R Function I create to be used within the
Kalman filter, and below this R function is my method for calling it within
the C++ code. To complicate matters further, the R Function calls a C++
Function. To clarify: The Kalman Filter C++ code calls an R function, and
this R Function calls an additional separate C++ function. (Code included
below)
As I iterate the Kalman filter it runs perfectly for anywhere from 30
minutes to six hours (and produces correct output when matched to the R
code), but it inevitably crashes. From past reading I have done, Dirk has
before mentioned that calling an R function many times within C++ is
usually not a good idea. As a result I suspect this is my issue (the error
codes vary, sometimes mentioning numerical errors, other times recursive
errors, other times random RCPP error codes -- I can document and provide
them if needed)
My biggest impasse is I cannot figure out a way to complete this without
calling the R optimize() function, and cannot find any RCPP optimize
functions to use instead. Thank you for reading.
*R Function (**IRNPV.CPP is the C++ function, which R optimizes the rate
parameter over):*
optimcpp<-function(CoupDate=1,coupNo=1,coup=1,price=1,rate=1)
{
m<-optimize(function(CoupDate,coupNo,coup,price,rate) *IRNPV.CPP*
(CoupDate=CoupDate,coupNo=coupNo,coup=coup,*rate*
)-price)^2,c(-0.05,0.2),tol=1e-20,CoupDate=CoupDate,coupNo=coupNo,coup=coup,price=price)
m$minimum
}
*Accessing the R environment within the C++ code:*
CPP.SRC <- '
using namespace arma;
Rcpp::Environment global = Rcpp::Environment::global_env();
Function optimizecpp = global["optimcpp"];
//Various matrix computations to calculate CD, CN, Co, & Pricex[0]
optimvec0 = optimizecpp(CD,CN,Co,pricex[0],Ra);
'
*IRNPV.CPP Function *(What the R Optimize() function optimizes 'rate' over
-- *Very* *Likely Unnecessary for Purposes of this Question)*
IR.NPV.TIPS.CBF.SRC <- '
using namespace arma;
double CD = Rcpp::as<double>(CoupDate);
double CN = Rcpp::as<double>(coupNo);
double RN = Rcpp::as<double>(rate);
double Co = Rcpp::as<double>(coup);
Rcpp::NumericVector LM;
mat DiscountFunc;
double length;
double price;
if (CN > 1) {
length = floor((((CD+0.5*(CN-1))-CD)/0.50))+1;
for (double i = CD; i <= (CD+0.5*(CN-1))+.05; i += 0.5) {
LM.insert(LM.end(),i);
}
DiscountFunc.set_size(LM.size(), 1);
DiscountFunc.fill(0);
double k = 0;
mat::row_iterator q = DiscountFunc.begin_row(0);
mat::row_iterator w = DiscountFunc.end_row((LM.size()-1));
for(mat::row_iterator i=q; i!=w; ++i)
{
(*i) = exp(-RN*LM[k]);
k = k + 1;
}
price = CD*Co*DiscountFunc[0];
for (int i=1; i<(LM.size()); ++i) {
price = price+0.5*Co*DiscountFunc[i];
}
}
else {
double DiscountFunc;
DiscountFunc = exp(-RN*CD);
price = (1+CD*Co)*DiscountFunc;
}
return Rcpp::wrap(price);
'
IRNPV.CPP <- cxxfunction(signature(CoupDate="NumericVector",
coupNo="NumericVector", coup="NumericVector",rate="NumericVector"),
IR.NPV.TIPS.CBF.SRC, include=CMATH, plugin="RcppArmadillo")
Thank you,
Simon
--
Simon Alexander Riddell
Economic Research RA
Federal Reserve Bank
_______________________________________________
Rcpp-devel mailing list
Rcpp-devel@lists.r-forge.r-project.org
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
