期权价值分析与风险分析的C++实现

risk_free_rate采用一年期shibor利率, expiration_date以距离到期日的自然日/365为准. 经过与东方财富网上的结果对比, 误差在0.1或者1%以内.

#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <sys/time.h>
using namespace std;

double sigma_call, sigma_put;
double value_call, value_put;
double price_call, price_put;
double delta_call, delta_put;
double gamma_call, gamma_put;
double vega_call,  vega_put;
double theta_call, theta_put;
double rho_call,   rho_put;

long getCurMicroSecond()
{
    struct timeval tm;
    gettimeofday(&tm, NULL);
    return ( (tm.tv_sec % 86400 + 28800) % 86400 )*1000000 + tm.tv_usec;
}


double normalDistribution(double d)
{
    return exp(-d*d/2)/sqrt(2*3.1415926);
}

double cumNormalDistribution(double d)
{
    if (d > 6)  return 1;
    if (d < -6) return 0;

    double gam = 0.231641900;
    double a1 = 0.319381530;
    double a2 = -0.356563782;
    double a3 = 1.781477937;
    double a4 = -1.821255978;
    double a5 = 1.330274429;

    double k = 1.0 / (1 + fabs(d) * gam);
    double n = k * (a1 + k * (a2 + k * (a3 + k * (a4 + k * a5))));
    n = 1 - normalDistribution(d) * n;
    if (d < 0) return 1.0 - n;

    return n;
}

double optionPricingCall(double s, double u, double e, double h, double r)
{
    double d1 = (log(u/s) + (r + (h*h)/2.0)*e)/(h*sqrt(e));
    double d2 = d1 - h*sqrt(e);
    return u*cumNormalDistribution(d1) - s*exp(-r*e) * cumNormalDistribution(d2);
}

double optionPricingPut(double s, double u, double e, double h, double r)
{
    double d1 = (log(u/s) + (r + (h*h)/2.0)*e)/(h*sqrt(e));
    double d2 = d1 - h*sqrt(e);
    return s*exp(-r*e) * cumNormalDistribution(-d2) - u*cumNormalDistribution(-d1);
}

double impliedVolatilityCall(double p,  double s, double u, double h, double e, double r)
{
    double sigma_pre = h;
    double price = optionPricingCall(s, u, e, sigma_pre, r);
    double sigma_next = sigma_pre;
    double sigma = sigma_pre;

    int count = 0;
    double step = 0.05;
    int MAX_CNT = 50;
    double ACCURACY = 1E-5; 

    if (price < p) {
        while (price < p && count < MAX_CNT) {
            sigma_pre = sigma_next;
            sigma_next = sigma_next + step;
            price = optionPricingCall(s, u, e, sigma_next, r);
            ++count;
        }
    }
    else if (price > p) {
        while (price > p && count < MAX_CNT) {
            sigma_next = sigma_pre;
            sigma_pre = sigma_pre - step;
            if (sigma_pre <= 0.0) {
                sigma_pre = 0.0;
                break;
            }
            price = optionPricingCall(s, u, e, sigma_next, r);
            ++count;
        }
    }
    count = 0;
    double diff = ((p-price) < 0 )? (price-p): (p-price);
    while ( (diff > ACCURACY) && (count < MAX_CNT)) {
        sigma = (sigma_pre + sigma_next) / 2;
        price = optionPricingCall(s, u, e, sigma, r);
        if (price < p) {
            sigma_pre = sigma;
        }
        else if (price > p) {
            sigma_next = sigma;
        }
        diff = ((p-price) < 0 )? (price-p): (p-price);
        ++count;
    }
    return sigma;
}

double impliedVolatilityPut(double opt_px,  double s, double u, double h, double e, double r)
{
    double sigma_pre = h;
    double price = optionPricingPut(s, u, e, sigma_pre, r);
    double sigma_next = sigma_pre;
    double sigma = sigma_pre;

    int count = 0;
    double step = 0.05;
    int MAX_CNT = 50;
    double ACCURACY = 1E-5; 

    if (price < opt_px) {
        while (price < opt_px && count < MAX_CNT) {
            sigma_pre = sigma_next;
            sigma_next = sigma_next + step;
            price = optionPricingPut(s, u, e, sigma_next, r);
            ++count;
        }
    }
    else if (price > opt_px) {
        while (price > opt_px && count < MAX_CNT) {
            sigma_next = sigma_pre;
            sigma_pre = sigma_pre - step;
            if (sigma_pre <= 0.0) {
                sigma_pre = 0.0;
                break;
            }
            price = optionPricingPut(s, u, e, sigma_next, r);
            ++count;
        }
    }
    count = 0;
    double diff = ((opt_px-price) < 0 )? (price-opt_px): (opt_px-price);
    while ( (diff > ACCURACY) && (count < MAX_CNT)) {
        sigma = (sigma_pre + sigma_next) / 2;
        price = optionPricingPut(s, u, e, sigma, r);
        if (price < opt_px) {
            sigma_pre = sigma;
        }
        else if (price > opt_px) {
            sigma_next = sigma;
        }
        diff = ((opt_px-price) < 0 )? (price-opt_px): (opt_px-price);
        ++count;
    }
    return sigma;
}

double Normal(double zz)
{ 
    //cdf of 0 is 0.5
    if (zz == 0) return 0.5;
    double z = zz;
    if (zz < 0) z = -zz;
    
    double p = 0.2316419;  
    double b1 = 0.31938153;
    double b2 = -0.356563782;
    double b3 = 1.781477937;
    double b4 = -1.821255978;
    double b5 = 1.330274428;
    
    double f = 1 / sqrt(2 * M_PI);
    double ff = exp(-pow(z, 2) / 2) * f;
    double s1 = b1 / (1 + p * z);
    double s2 = b2 / pow((1 + p * z), 2);
    double s3 = b3 / pow((1 + p * z), 3);
    double s4 = b4 / pow((1 + p * z), 4);
    double s5 = b5 / pow((1 + p * z), 5);
    
    double  sz = ff * (s1 + s2 + s3 + s4 + s5); 
    double rz; 
    if (zz < 0) rz = sz;
    if (zz > 0) rz = (1 - sz);
    return rz;
}

double valueCall(double s, double u, double sigma, double r, double e)
{ 
     double ls = log(u);
     double lx = log(s);
     double t = e;
     double sigma2 = pow(sigma, 2);
     double n = (ls - lx + r * t + sigma2 * t / 2);
     double sqrtT = sqrt(e);
     double d = sigma * sqrtT;
     double d1 = n / d;
     double d2 = d1 - sigma * sqrtT;
     double nd1 = Normal(d1);
     double nd2 = Normal(d2);
     return u * nd1 - s * exp(-r * t) * nd2;
}
    
double valuePut(double s, double u, double sigma, double r, double e)
{
     double ls = log(u);
     double lx = log(s);
     double t = e;
     double sigma2 = pow(sigma, 2);
     double n = (ls - lx + r * t + sigma2 * t / 2);
     double sqrtT = sqrt(e);
     double d = sigma * sqrtT;
     double d1 = n / d;
     double d2 = d1 - sigma * sqrtT;
     double nd1 = Normal(d1);
     double nd2 = Normal(d2);
     return s * exp(-r * t) * (1 - nd2) - u * (1 - nd1);
}

void greeksCall(double s, double u, double e, double sigmas, double r)
{
    double d1 = (log(u / s)  + ((r + sigmas * sigmas * e) / 2.0)) / (sigmas * sqrt(e));
    double d2 = d1 - sigmas * sqrt(e);
    double c1 = u * cumNormalDistribution (d1);
    delta_call = exp(-r * e) * cumNormalDistribution(d1);
    vega_call = u * sqrt(e) * normalDistribution(d1) * exp(-r * e);
    gamma_call = exp(-r * e) * normalDistribution(d1) / (u * sigmas * sqrt(e));
    theta_call = exp(-r * e) * (-u * normalDistribution(d1) * sigmas / (2 * sqrt(e)) 
                    - r * s * cumNormalDistribution(d2) + r * u * cumNormalDistribution(d1));
    rho_call = e * s * exp (-r * e) * cumNormalDistribution (d2);

    price_call = u * cumNormalDistribution(d1) - s * exp(-r * e) * cumNormalDistribution(d2);
}


void greeksPut(double s, double u, double e, double sigmas, double r)
{
    double d1 = (log(u / s)  + ((r + sigmas * sigmas * e) / 2.0)) / (sigmas * sqrt(e));
    double d2 = d1 - sigmas * sqrt(e);
    double c1 = u * cumNormalDistribution (d1);
    delta_put = exp(-r * e) * (cumNormalDistribution(d1) - 1);
    vega_put = u * sqrt(e) * normalDistribution(d1) * exp(-r * e);
    gamma_put = exp(-r * e) * normalDistribution(d1) / (u * sigmas * sqrt(e));
    theta_put = exp(-r * e) * (-u * normalDistribution(d1) * sigmas / (2 * sqrt(e)) 
                    - r * s * cumNormalDistribution(d2) + r * u * cumNormalDistribution(d1));
    rho_put = e * s * exp (-r * e) * cumNormalDistribution (d2);

    price_put = u * cumNormalDistribution(d1) - s * exp(-r * e) * cumNormalDistribution(d2);
}

int main(){
    double p = 0.3591;//option_price
    double u = 3.080;//underlying_price
    double s = 2.75;//strike
    double r = 0.03032;//risk_free_rate
    double e = 0.197260;//expiration_date
    double h = 0;//history_volatility

    sigma_call = impliedVolatilityCall(p, s, u, h, e, r);
    sigma_put  = impliedVolatilityPut(p, s, u, h, e, r);

    value_call = valueCall(s, u, sigma_call, r, e);
    value_put  = valuePut(s, u, sigma_put, r, e);

    greeksCall(s, u, e, sigma_call, r);
    greeksPut(s, u, e, sigma_put, r);

    printf("option_price=%lf\nunderlying_price=%lf\nstrike=%lf\nrisk_free_rate=%lf\nexpiration_date=%lf\nhistory_volatility=%lf\n"
        , p, u, s, r, e, h);
    printf("sigma_call=%lf\ndelta_call=%lf\ngamma_call=%lf\nvega_call=%lf\ntheta_call=%lf\nrho_call=%lf\nprice_call=%lf\nvalue_call:%lf\n"
        , sigma_call, delta_call, gamma_call, vega_call, theta_call, rho_call, price_call, value_call);
    printf("sigma_put=%lf\ndelta_put=%lf\ngamma_put=%lf\nvega_put=%lf\ntheta_put=%lf\nrho_put=%lf\nprice_put=%lf\nvalue_put:%lf\n"
        , sigma_put, delta_put, gamma_put, vega_put, theta_put, rho_put, price_put, value_put);
    return 0;
}