-1

ここで助けが必要です。コードの複雑さを許してください。基本的に、二分法を使用して値「Theta」を見つけ、i がインクリメントするたびに探しています。シータがわかっていればすべての計算が正常に機能することがわかっており、コードを実行してすべての値を単純に計算するだけですが、while ループと二分法を導入してコードにシータを近似させると、うまくいきません。正しく実行されるようです。while ループの設定が間違っていると思います....

#include <math.h>
#include <iostream>
#include <vector>
#include <iomanip>
#include <algorithm>    // std::max

using namespace std;

double FuncM(double theta, double r, double F, double G, double Gprime, double d_t, double sig);

double FuncM(double theta, double r, double F, double G, double Gprime, double d_t, double sig)
{
    double eps = 0.0001;
    return ((log(max((r + (theta + F - 0.5 * G * Gprime ) * d_t), eps))) / sig); 
}

double FuncJSTAR(double m, double x_0, double d_x);

double FuncJSTAR(double m, double x_0, double d_x)
{
    return (int(((m - x_0) / d_x)+ 0.5));
}

double FuncCN(double m, double x_0, double j, double d_x);

double FuncCN(double m, double x_0, double j, double d_x)
{
    return (m - x_0 - j * d_x);
}

double FuncPup(double d_t, double cn, double d_x);

double FuncPup(double d_t, double cn, double d_x)
{
    return (((d_t + pow(cn, 2.0)) / (2.0 * pow(d_x, 2.0))) + (cn / (2.0 * d_x)));
}

double FuncPdn(double d_t, double cn, double d_x);

double FuncPdn(double d_t, double cn, double d_x)
{
    return (((d_t + pow(cn, 2.0)) / (2.0 * pow(d_x, 2.0))) - (cn / (2.0 * d_x)));
}

double FuncPmd(double pd, double pu);

double FuncPmd(double pd, double pu)
{
    return (1 - pu - pd);
}

int main()
{
    const int Maturities = 5;
    const double EPS = 0.00001;

    double TermStructure[Maturities][2] = {
        {0.5 , 0.05},
        {1.0 , 0.06},
        {1.5 , 0.07},
        {2.0 , 0.075},
        {3.0 , 0.085}               };

//--------------------------------------------------------------------------------------------------------
    vector<double> Price(Maturities);

    double Initial_Price = 1.00;

    for (int i = 0; i < Maturities; i++)
        {
            Price[i] = Initial_Price * exp(-TermStructure[i][1] * TermStructure[i][0]);
        }
//--------------------------------------------------------------------------------------------------------
    int j_max = 8;
    int j_range = ((j_max * 2) + 1);
//--------------------------------------------------------------------------------------------------------
// Set up vector of possible j values
    vector<int> j_value(j_range);

    for (int j = 0; j < j_range; j++)
        {
            j_value[j] = j_max - j;
        }
//--------------------------------------------------------------------------------------------------------
    double dt = 0.5;
    double dx = sqrt(3 * dt);
    double sigma = 0.15;
    double mean_reversion = 0.2; // "a" value
//--------------------------------------------------------------------------------------------------------
    double r0 = TermStructure[0][1]; // Initialise r(0) in case no corresponding dt rate in term structure
//--------------------------------------------------------------------------------------------------------
    double x0 = log(r0) / sigma;
//--------------------------------------------------------------------------------------------------------
    vector<double> r_j(j_range); // rate at each j
    vector<double> F_r(j_range);
    vector<double> G_r(j_range);
    vector<double> G_prime_r(j_range);

    for(int j = 0; j < j_range; j++)
        {
            if (j == j_max)
                {
                    r_j[j] = r0;
                }
            else
                {
                    r_j[j] = exp((x0 + j_value[j]*dx) * sigma);
                }

            F_r[j] = -mean_reversion * r_j[j];
            G_r[j] = sigma * r_j[j];
            G_prime_r[j] = sigma;
        }
//--------------------------------------------------------------------------------------------------------

    vector<vector<double>> m((j_range), vector<double>(Maturities));
    vector<vector<int>> j_star((j_range), vector<int>(Maturities));
    vector<vector<double>> Central_Node((j_range), vector<double>(Maturities));
    vector<double> Theta(Maturities - 1);

    vector<vector<double>> Pu((j_range), vector<double>(Maturities));
    vector<vector<double>> Pd((j_range), vector<double>(Maturities));
    vector<vector<double>> Pm((j_range), vector<double>(Maturities));

    vector<vector<double>> Q((j_range), vector<double>(Maturities));// = {}; // Arrow Debreu Price. Initialised all array values to 0
    vector<double> Q_dt_sum(Maturities);// = {}; // Sum of Arrow Debreu Price at each time step. Initialised all array values to 0

//--------------------------------------------------------------------------------------------------------

    double Theta_A, Theta_B, Theta_C;
    int JSTART;
    int JEND;
    int TempStart;
    int TempEnd;
    int max;
    int min;
    vector<vector<int>> Up((j_range), vector<int>(Maturities));
    vector<vector<int>> Down((j_range), vector<int>(Maturities));

//  Theta[0] = 0.0498039349327417;
//  Theta[1] = 0.0538710670441647;
//  Theta[2] = 0.0181648634139392;
//  Theta[3] = 0.0381183886467521;

    for(int i = 0; i < (Maturities-1); i++)
        {
            Theta_A = 0.00;
            Theta_B = TermStructure[i][1];
            Q_dt_sum[0] = Initial_Price;
            Q_dt_sum[i+1] = 0.0;

            while (fabs(Theta_A - Theta_B) >= 0.0000001)
                {
                max = 1;
                min = 10;

                if (i == 0)
                    {
                        JSTART = j_max;
                        JEND = j_max;
                    }
                else
                    {
                        JSTART = TempStart;
                        JEND = TempEnd;
                    }   

                for(int j = JSTART; j >= JEND; j--)
                    {

                        Theta_C = (Theta_A + Theta_B) / 2.0; // If Theta C is too low, the associated Price will be higher than Price from initial term structure. (ie P(Theta C) > P(i+2) for Theta C < Theta)
                        // If P_C > P(i+2), set Theta_B = Theta_C, else if P_C < P(i+2), set Theta_A = Theta_C, Else if P_C = P(i+2), Theta_C = Theta[i]
                        //cout << Theta_A << "  " << Theta_B << "       " << Theta_C << endl;
                        m[j][i] = FuncM(Theta[i], r_j[j], F_r[j], G_r[j], G_prime_r[j], dt, sigma);
                        j_star[j][i] = FuncJSTAR(m[j][i], x0, dx);
                        Central_Node[j][i] = FuncCN(m[j][i], x0, j_star[j][i], dx);
                        Pu[j][i] = FuncPup(dt, Central_Node[j][i], dx);
                        Pd[j][i] = FuncPdn(dt, Central_Node[j][i], dx);
                        Pm[j][i] = FuncPmd(Pd[j][i], Pu[j][i]);

                        for (int p = 0; p < j_range; p++)
                            {
                                Q[p][i] = 0; // Clear Q array
                            }

                        Q[j_max][0] = Initial_Price;
                        Q[j_max -(j_star[j][i]+1)][i+1]     = Q[j_max - (j_star[j][i]+1)][i+1] + Q[j][i] * Pu[j][i] * exp(-r_j[j] * dt); 
                        Q[j_max -(j_star[j][i]  )][i+1]     = Q[j_max - (j_star[j][i]  )][i+1] + Q[j][i] * Pm[j][i] * exp(-r_j[j] * dt); 
                        Q[j_max -(j_star[j][i]-1)][i+1]     = Q[j_max - (j_star[j][i]-1)][i+1] + Q[j][i] * Pd[j][i] * exp(-r_j[j] * dt);
                    }

                    for (int j = 0; j < j_range; j++)
                        {
                            Up[j][i] = j_star[j][i] + 1;
                            Down[j][i] = j_star[j][i] - 1;

                            if (Up[j][i] > max)
                                {
                                    max = Up[j][i];
                                }

                            if ((Down[j][i] < min) && (Down[j][i] > 0))
                                {
                                    min = Down[j][i];
                                }
                        }

                        TempEnd = j_max - (max);
                        TempStart = j_max - (min);

                        for (int j = 0; j < j_range; j++)
                            {
                                Q_dt_sum[i+1] = Q_dt_sum[i+1] + Q[j][i] * exp(-r_j[j] * dt);                
                                cout << Q_dt_sum[i+1] << endl;
                            }


                        if (Q_dt_sum[i+1] == Price[i+2])
                            {
                                Theta[i] = Theta_C;
                                break;
                            }

                        if (Q_dt_sum[i+1] > Price[i+2])
                            {
                                Theta_B = Theta_C;
                            }
                        else if (Q_dt_sum[i+1] < Price[i+2])
                            {
                                Theta_A = Theta_C;
                            }

                }
        cout << Theta[i] << endl;       
        }
    return 0;
}
4

1 に答える 1

0

わかりました、悪いです。間違って呼び出された値がありました。すべて良い。

于 2016-05-05T04:43:59.843 に答える