#--------------------------#
#        Chapter 1:        #
#   Motivating Example 2   #
#--------------------------#


mu = 2
lam1 = -1
lam2 = +1/2
lam3 = +3/2

beta = 200                             
dt = 0.1
tk = seq(0,beta,by=dt)
nt = length(tk)

w1 = exp(tk*lam1^2/2)
w2 = exp(tk*lam2^2/2)
w3 = exp(tk*lam3^2/2)
EGth = ( w1*cos(mu*lam1*sqrt(tk)) + w2*cos(mu*lam2*sqrt(tk)) + w3*cos(mu*lam3*sqrt(tk)) )/(w1+w2+w3)


set.seed(123456)
N = 100000                                   # number MC simus


#  untransformed plain Monte Carlo:

x = rnorm(N)
xsqrbeta = outer(x,sqrt(tk))
Z = exp(lam1*xsqrbeta) + exp(lam2*xsqrbeta) + exp(lam3*xsqrbeta)
G = exp(mu^2/2) * cos(mu*x) 
G = outer(G,rep(1,nt))
EG0 = colSums(G*Z)/colSums(Z)


# picture:

info = "untransformed vs. Girsanov transformed Monte Carlo, 
green = exact, blue = transformed MC, red = untransformed MC" 

plot(tk, EG0, col="red", cex=0.4, main=info, font.main=10, xlab="beta", ylab="< G >" )
points(tk, EGth, col="green", cex=1.0 )



#  Girsanov transformed Monte Carlo (may take 1 or 2 minutes):

Z1 = function( x )
{
  e1 = exp(lam1*x)
  e2 = exp(lam2*x)
  e3 = exp(lam3*x)
  res = (lam1*e1 + lam2*e2 + lam3*e3)/(e1+e2+e3)
  return(res)
}
Z2 = function( x )
{
  e1 = exp(lam1*x)
  e2 = exp(lam2*x)
  e3 = exp(lam3*x)
  res = (lam1^2*e1 + lam2^2*e2 + lam3^2*e3)/(e1+e2+e3)
  return(res)
}

x = matrix(0, nrow=N, ncol=nt)
intZ2 = matrix(0, nrow=N, ncol=nt)
prob = matrix(0, nrow=N, ncol=nt)
for(k in 2:nt)
{
  phi = rnorm(N)
  x[,k] = x[,k-1] + Z1(x[,k-1])*dt + sqrt(dt)*phi
  intZ2[,k] = intZ2[,k-1] + 1/2*Z2(x[,k-1])*dt
  maxintZ2 = max( intZ2[,k] )
  prob[,k] = exp( intZ2[,k] - maxintZ2 )
} 
G = exp(mu^2/2) * cos(mu*x/outer(rep(1,N),sqrt(tk))) 
EG = colSums(G*prob)/colSums(prob)
points(tk, EG, cex=0.4, col="blue" )