GSL과의 C - (수학적) 통합은 세그먼트 화 오류를 야기합니다.

Question

그래서 저는 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); 

} 

Answer 1

귀하의 hestonIntegrand1 (나는 hestonIntegrand2 생각) 함수는 첫 번째 매개 변수로 복잡한 번호를 가지고 있지만, gsl_function 만 double받는 함수를 할 수 있습니다. 가장 높은 경고 수준을 설정 한 경우 호환되지 않는 포인터 할당을 확인해야합니다.

복잡한 함수의 각 구성 요소에 대해 별도의 함수를 만들고 실제 구성 요소에 대해 한 번, 허수 구성 요소에 대해 한 번 두 번 통합하여이 제한을 해결할 수 있습니다. 결과를 함께 결합하면 복잡한 결과를 얻을 수 있습니다.