lambda<-0.0375
sigma<-0.017607
sigmaJ<-0.097368
A<-0.0234
xbar<-3.71903
for(i in 1:10000){
    sim[i,]<-GenerateOnePath(100,ini,1,lambda,sigma,sigmaJ,A,xbar)
    }

/**** Simulating one path for given parameters ****/

GenerateOnePath<-function(NPE,ini,year,lambda,sigma,sigmaJ,A,xbar){
  path<-rep(0,NPE)
  path[1]<-ini    
  slope<-exp(-A*year*365/NPE)
  delta<-year*365/NPE
  intercept<-xbar*(1-slope)
  csigma<-sigma*sqrt((1-exp(-2*A*year*365/NPE))/2/A)
  for(i in 2:NPE){
                  NOJ<-rpois(1,lambda*delta)
                  jump<-0
                  if(NOJ>0){    
                            for(j in 1:NOJ){
                        jump<-jump+exp(-A*delta*(1-runif(1,min=0,max=1)))*sigmaJ*rnorm(1)
                                           }
                           }
                  path[i]<-slope*path[i-1]+intercept+csigma*rnorm(1)+jump
                  }
path
}

GenerateOnePath1<-function(NPE,ini,year,sigma,){
  path<-rep(0,NPE)
  path[1]<-ini    
  delta<-year*365/NPE
  for(i in 2:NPE){
                  path[i]<-path[i-1]+(sigma*sqrt(delta)*rnorm(1)
                  }
path
}

/****** Decompose continous state space into discete states for optimization   
        We will denote finer resolutions for states around long-term mean, and 
        coarser resolution for states far away from long-term mean             ******/

/*****  sd         ---- stands for standard deviation of the stochastic process
        lmean      ---- stands for long-term mean of the stochstic process
        NumofStates---- stands for number of states intended to decompose,we assume
                        that the number is divisble by 4                       ******/


GetStates<-function(sd,lmean,NumofStates){
      
               result<-rep(0,NumofStates-1)
               n<-NumofStates/4
               result[1]<-lmean-14*sd;               
               for(i in 2:n) {
                      result[i]<-result[1]+(i-1)*32*sd/NumofStates
                             }
               for(i in 1:n) {
                      result[i+n]<-result[n]+i*24*sd/NumofStates
                             }
               for(i in 1:n) {
                      result[i+2*n]<-result[2*n]+i*24*sd/NumofStates
                             }
               m<-n-1   
               for(i in 1:m) {
                      result[i+3*n]<-result[3*n]+i*32*sd/NumofStates
                             }
               result
}



/***** Definition of h function, based on strike price
       Strike   ----- strike price defined on option 
        x       ----- underlying state                  ******/

h<-function(x,Strike){
      if (exp(x)>=Strike) value<-exp(x)-Strike
      else value<-0
}


/*******
       States     ------  State-space limits, a vector, with length 40 right now
       TwoStages  ------  Collection of two consecutive points for all path at time
                          t=2,...,NPE-1.
       Iterations ------  Numbers of paths simulated. 
       Strike     ------  strike price defined on option                              *******/



GetFirstV<-function(Iterations,V,States){
       
     NumofStates<-length(States)+1
     SumofV<-rep(0,NumofStates)
     FirstV<-rep(0,NumofStates)
     Frequency<-rep(0,NumofStates)   
     for(i in 1:Iterations){

                 j<-sum(States<V)+1
                 Frequency[j]<- Frequency[j]+1          
                 SumofV[j]<-SumofV[j]+V[i]

                           }
     for(j in 1:NumofStates){
               if(Frequency[j]>0){
                    FirstV[j]<-SumofV[j]/Frequency[j]
                                 }
               else FirstV[j]<-0
                            }
     FirstV
}
/**** Get phat[j,k] and hhat[j] for every time t=1,...,NPE   ****/


GetPhatAndHhat<-function(Iterations,TwoStages,States,Strike){
       
     NumofStates<-length(States)+1
     Frequency<-array(0,dim=c(NumofStates,NumofStates)) 
     phat<-array(0,dim=c(NumofStates,NumofStates)) 
     Sumofh<-rep(0,NumofStates)
     hhat<-rep(0,NumofStates)
   
     for(i in 1:Iterations){

                 j<-sum(States<TwoStages[i,1])+1
                 k<-sum(States<TwoStages[i,2])+1
                 Frequency[j,k]<- Frequency[j,k]+1          
                 Sumofh[j]<-Sumofh[j]+h(TwoStages[i,1],Strike)

                           }

     for(j in 1:NumofStates){

                 for(k in 1:NumofStates){
                         
                        if(sum(Frequency[j,])>0){              

                           phat[j,k]<-Frequency[j,k]/sum(Frequency[j,])
                                         
                                                 }
             
                        else phat[j,k]<-0 
                                    
                                        }                 
                            
                            }

     for(j in 1:NumofStates){

                 if(sum(Frequency[j,])>0){              

                 hhat[j]<-Sumofh[j]/sum(Frequency[j,])
                                         
                                         }
                 else hhat[j]<-0                  
                            
                            }
     result<-list(phat,hhat)
     result
 }

/***** One stage maxmization given next time maxmization values V,phat and hhat *****/

/*****      V        ------ stands for vector Vhat[i+1]                                         
            Phat     ------ stands for matrix phat[i]
            hhat     ------ stands for vector hhat[i]  
            discount ------ stands for discount factor for present value        *****/

GetMax<-function(V,phat,hhat,discount){
         
         NumofStates<-length(hhat)
         result<-rep(0,NumofStates) 
         Temp<-phat%*%V
         
         for(i in 1:NumofStates){

                    result[i]<-max(hhat[i],discount*Temp[i,1])
 
                                }
         result
}

Pricing<-function(sim,Iterations,NPE,year,sd,lmean,NumofStates,Strike,interest){
          States<-GetStates(sd,lmean,NumofStates)
          discount<-exp(-interest*year/NPE)
          V<-h(sim[,NPE],Strike)
          FirstV<-GetFirstV(Iterations,V,States)
          for(i in 1:99){
          firstcol<-NPE-i
          nextcol<-firstcol+1
          hatlist<-GetPhatAndHhat(Iterations,sim[,firstcol:nextcol],States,Strike)                  
          FirstV<-GetMax(FirstV,hatlist[[1]],hatlist[[2]],discount)
                        }
          FirstV
}

Num<-100
sim1<-array(0,dim=c(Num,10))
for(i in 1:Num) {
     sim1[i,]<-GenerateOnePath1(10,ini,1,0.025)
}

price<-Pricing(sim1,100,10,1,0.025,ini,16,44,0.0246)