C++ whileループは二分法を使用しています。休憩に関するヘルプ

Question

私はここでいくつかの助けが必要です。コードの複雑さを許してください。基本的には、私は二分法を使って値 "シータ"とそれぞれのiの増分を探しています。 私はThetaを知っていればすべての計算がうまく動作することを知っています。コードを実行してすべての値を計算するだけですが、Wh​​ileループと二分法を導入してThetaを近似すると、それが正しく動作するように見える。私は私の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; 
} 

Answer 1

Ok、私の悪いです。私は値が間違って呼ばれていた。 良い