そこで、Cでオプションモデルの統合を行おうとしています。関数はMatlabから呼び出され、mexファイルとしてコンパイルされます。関数hestonIntegrand1は、gsl_integration_qaiguの外部で呼び出されたときに値を返すため、この統合関数の内部でエラーが発生している必要があります。何が問題なのですか？それは単なる数学の問題ですか？
#include "mex.h"
#include "math.h"
#include "complex.h"
#include <gsl/gsl_integration.h>
#define MAX_WORKSPACE_SIZE 1000
struct hestonIntegrandParams
{
double _Complex K;
//double X;
double St;
double tau;
double theta;
double kappa;
double sigmaV;
double sigma;
double rho;
double gamma;
double r;
};
double hestonIntegrand1 (double _Complex phi, void * p) {
/*
double _Complex im = 0.0f + 1.0f * _Complex_I ;
double _Complex re = 1.0f + 0.0f * _Complex_I ;
*/
struct hestonIntegrandParams * params = (struct hestonIntegrandParams *)p;
double _Complex K = (params->K);
//double X = (params->X);
double St = (params->St);
double tau = (params->tau);
double theta = (params->theta);
double kappa = (params->kappa);
double sigmaV = (params->sigmaV);
double sigma = (params->sigma);
double rho = (params->rho);
double gamma = (params->gamma);
double r = (params->r);
double x = log (St);
double u1 = 0.5f;
//double b1 = kappa+lambda-rho*sigma;
// under EMQ lambda falls out
double b1 = kappa-rho*sigma;
double _Complex i = 0.0f + 1.0f * _Complex_I;
double _Complex d1 = sqrt(pow((rho*sigma*phi*i-b1),2) - (sigma*sigma)*(2*u1*phi*i-phi*phi) );
double _Complex g1 = (b1-rho*sigma*phi*i+d1) / (b1-rho*sigma*phi*i-d1);
double _Complex D1 = (b1-rho*sigma*phi*i+d1)/(sigma*sigma) * ((1.0f-exp(d1*tau))/(1.0f-g1*exp(d1*tau)));
double _Complex C1 = r*phi*i*tau + (kappa*theta/(sigma*sigma))* (
(b1-rho*sigma*phi*i+d1)*theta -
2 - log( (1- g1*exp(d1*theta) ) / (1-g1) )
);
double _Complex f1 = exp(C1 + D1*sigmaV + i*phi*x);
double _Complex Re = 1.0f + 0.0f * _Complex_I;
double _Complex integrand = Re*(exp(-i*phi*log(K))*f1)/(i*phi);
// returning just the real part
return (creal(integrand));
}
// -callHestoncf(u, X, v, r, q, v0, vT, rho, k, sigma, implVol )
double hestonCallOption ( double St, double K,
double tau, double r, double theta, double kappa, double sigma, double sigmaV,
double rho, double gamma) {
gsl_integration_workspace *work_ptr =
gsl_integration_workspace_alloc (MAX_WORKSPACE_SIZE);
gsl_integration_workspace *work_ptr2 =
gsl_integration_workspace_alloc (MAX_WORKSPACE_SIZE);
double lower_limit = 0.0f; /* start integral from 1 (to infinity) */
double abs_error = 1.0e-8; /* to avoid round-off problems */
double rel_error = 1.0e-8; /* the result will usually be much better */
/* the result from the integration
and the estimated errors from the integration
*/
double result1, result2, error1, error2;
double alpha = 2.0;
double expected = -0.16442310483055015762; /* known answer */
gsl_function F1;
gsl_function F2;
struct hestonIntegrandParams params = {
K, St, tau, theta, kappa, sigmaV, sigma, rho, gamma, r
};
F1.function = &hestonIntegrand1;
F1.params = ¶ms;
F2.function = &hestonIntegrand2;
F2.params = ¶ms;
mexWarnMsgTxt("YOLO");
gsl_integration_qagiu (&F1, lower_limit,
abs_error, rel_error, MAX_WORKSPACE_SIZE, work_ptr, &result1,
&error1);
gsl_integration_qagiu (&F2, lower_limit,
abs_error, rel_error, MAX_WORKSPACE_SIZE, work_ptr2, &result2,
&error2);
mexWarnMsgTxt("YOLO2");
gsl_integration_workspace_free (work_ptr);
gsl_integration_workspace_free (work_ptr2);
// df = discount factor
double df = exp(-r*tau);
//StP1 ? KP(t, T)P2
double price = St*result1-K*df*result2;
//double price = hestonIntegrand1(1, ¶ms);
//double price = 2.0f;
return price;
}
/* the gateway function */
void mexFunction( int nlhs, mxArray *plhs[],
int nrhs, const mxArray *prhs[])
{
double *y,*z;
double x;
mwSize mrows,ncols;
/* check for proper number of arguments */
/* NOTE: You do not need an else statement when using mexErrMsgTxt
within an if statement, because it will never get to the else
statement if mexErrMsgTxt is executed. (mexErrMsgTxt breaks you out of
the MEX-file) */
char buf[256];
snprintf(buf, sizeof buf, "11 inputs required, but only %i provided", nrhs);
if(nrhs!=1)
mexErrMsgTxt(buf);
if(nlhs!=1)
mexErrMsgTxt("1 output required.");
/* check to make sure the first input argument is a scalar */
/*
if( !mxIsDouble(prhs[0]) || mxIsComplex(prhs[0]) ||
mxGetN(prhs[0])*mxGetM(prhs[0])!=1 ) {
mexErrMsgTxt("Input x must be a scalar.");
}
*/
/* set the output pointer to the output matrix */
plhs[0] = mxCreateDoubleMatrix(mrows,ncols, mxREAL);
/* create a C pointer to a copy of the output matrix */
//z = mxGetPr(plhs[0]);
mxArray *xData;
double *xValues, *outArray;
//Copy input pointer x
xData = prhs[0];
xValues = mxGetPr(xData);
double price = hestonCallOption ( xValues[0] , xValues[1], xValues[2] ,
xValues[3], xValues[4], xValues[5], xValues[6], xValues[7],
xValues[8], xValues[9] );
//Allocate memory and assign output pointer
plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL); //mxReal is our data-type
outArray = mxGetPr(plhs[0]);
outArray[0] = price;
/* call the C subroutine */
//xtimesy(x,y,z,mrows,ncols);
}