//
// simu2.cc  Stochastic Volatility Option Pricing Model
//
#include <math.h>
#include <stdlib.h>
#include <iostream.h>
#include "simu2.h"
#define unifrand() ((double) rand()/((double) RAND_MAX))
#define MAX(a,b) (((a) > (b)) ? (a) : (b))

svprice::svprice(double ttm,int nper,double r,double initstd, double
    meanstd, double k, double sigstd) {
  npers = nper;
  tinc = ttm/(double) nper;
  r1per = 1.0 + r*tinc;
  sigma0 = initstd;
  meansigma = meanstd;
  sigmasigma = sigstd;
  kappa = k;
  c0 = kappa * tinc * meansigma;
  c1 = 1.0 - kappa * tinc;
  c2 = 1.0 - sigmasigma * sqrt(tinc)*0.5*sqrt((double) 12);
  c3 = sigmasigma * sqrt(tinc) * sqrt((double) 12);
  c4 = sqrt(tinc)*sqrt((double) 12);}

double svprice::eurcall(double stock,double strike,long int nsimu) {
  long int i,n;
  double x;
  x=0.0;
  for(n=0;n<nsimu;n++) {
    stockP = stock;
    sigma = sigma0;
    for(i=0;i<npers;i++) {
      stockP *= stocktotret();
      }
    x += MAX(stockP-strike,0);}
  return(x/((double) nsimu * pow(r1per,(double) npers)));}

double svprice::stocktotret() {
  double stockret;
//
//  The following straightforward computations are algebraically the same as
//  the ones used below but are much slower because many more calculations are
//  performed in each pass through the loop.
//
//  stockret = r1per + sigma * sqrt(tinc) * (unifrand()-0.5) * sqrt((double) 12);
//  sigma = (kappa*tinc * meansigma + (1.0 - kappa * tinc) * sigma) *
//    (1 + sigmasigma * sqrt(tinc) * (unifrand()-0.5) * sqrt((double) 12));
//  return(stockret);}
//
  stockret = r1per + sigma * c4 * (unifrand() - 0.5);
  sigma = (c0 + c1 * sigma) * (c2 + c3 * unifrand());
  return(stockret);}