1

所以我正在尝试在 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 = &params;
    F2.function = &hestonIntegrand2;
    F2.params = &params;

    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, &params);
    //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);

}
4

1 回答 1

1

您hestonIntegrand1（我相信hestonIntegrand2）函数将复数作为第一个参数，但gsl_function只允许使用double. 如果您启用了最高警告级别，您应该会看到不兼容的指针分配。

您可以通过为复杂函数的每个组件创建单独的函数并集成两次来解决此限制，一次用于实部，一次用于虚部。将结果组合在一起应该会给您带来复杂的结果。

于 2013-03-15T22:01:48.607 回答