/*@@@@@@@@@@@@@@@@@@@@@@@@ LocalMC.cpp @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
 Random Option Pricing
 Summer 2002
 Authors: A. Korobchevsky
          S. Ferrando
 Ryerson University ALL RIGHTS RESERVED
 
 Purpose: Discrete Binomial implementation of Local Monte Carlo Algorithm
 Note: The algorithm performs shift by any value in [1..N].
       When the shift parameter equals 1 the algorithm is equivalent to non-shift
       local Monte Carlo  
 @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/

#include <cmath>

#include "../Und/Binom.h" 
#include "LocalMC.h" 




/*@@@@@@@@@@@@@@@@@@@@@@@@  calcResults  @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
 Function: calcResults
 Purpose: Compute the option price and it's standard error
 Param: Underlying - underlying function object
        Derivative - derivative function object
	double& - Option value (output param.)  
	double& - Standard error (output param.)
 Precondition: Underlying, Derivative and RNG parameters initialized
 Note: Due to rounding errors caused by separation of simulations into classes 
       the total number of simulations might not equal to M. It must be computed by 
       using sum of simulations for each class.
 @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/
void LocalMC::calcResults(Underlying* binom, Derivative* ec, vector<double>& cv, vector<double>& se)
{ double sum_CT = 0;
  int sumSims = 0;  // Total number of path including the shifted ones

  double *sum_i = new double[binom->N + 1];  // Sum of payoffs in each class
  int* counter = new int[binom->N + 1];  // Number of paths in each class
  double* sum = new double[binom->N + 1];  // Sum of option values in each class
  double* error = new double[binom->N + 1];  // Standard error in each class

  for (int k = 0; k <= binom->N; k++)
  { sum_i[k] = sum[k] = error[k] = 0.0;
    counter[k] = 0;
  }

  for (int Mi = 0; Mi < M; Mi += gap) 
  { sum_CT = 0;

    for (int k = 0; k <= binom->N; k++)
    { int simulations = (int)((Mi + gap) * logProbability(1 - binom->P, binom->N, k)) - 
                          (int)(Mi * logProbability(1 - binom->P, binom->N, k));
      sumSims += simulations * shift;   
      for (int j = 1; j <= simulations; j++)
      { double* path = new double[binom->N];
        path = binom->generatePath(binom->N, k);

        double ST = S;
        ec->init();
        ec->update(ST); 
        for (int i = 0; i < binom->N; i++)
        { ST *= path[i];
          ec->update(ST); 
        }
        double CT = ec->payOff(ST);
        sum_i[k] += CT;
        if (counter[k] > 1)
          error[k] = ((counter[k] - 2) * error[k] + 
							(counter[k] - 1.0) / counter[k] * pow(ST - sum[k] / (counter[k] - 1), 2)) / (counter[k] - 1); 
		  counter[k]++;
		  sum[k] += ST;     

// BEGIN Shift
        for (int shifttimes = 0; shifttimes < shift - 1; shifttimes++)
        { ec->sum_St = (ec->sum_St - S) / path[shifttimes] + ST;
          if (ST > S)
            ec->min_St = MIN(ec->min_St * path[binom->N - shifttimes - 1], S);
          else
            ec->min_St = MIN(ec->min_St / path[shifttimes], ST);

          double CT = ec->payOff(ST);
          sum_i[k] += CT;  

          if (counter[k] > 1)
            error[k] = ((counter[k] - 2) * error[k] + 
							  (counter[k] - 1.0) / counter[k] * pow(ST - sum[k] / (counter[k] - 1), 2)) / (counter[k] - 1); 
		    counter[k]++;
		    sum[k] += ST;     
        }
// END Shift

      }  
    }

    if (Mi > M0)
    { for (int i = 0; i <= binom->N; i++)
        sum_CT += sum_i[i];
      cv.push_back(sum_CT/sumSims * exp(-binom->r * binom->T));

// BEGIN Error
      // Error computation
      double varTotal = 0;
		for (int k = 0; k <= binom->N; k++)
        if (counter[k] > 0)
          varTotal += pow(logProbability(1 - binom->P, binom->N, k), 2) * error[k] / counter[k];      
// END Error    
   
      se.push_back(exp(-binom->r * binom->T / binom->N) * sqrt(varTotal));
    }
  }
} 
