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


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

mu = 2
EGth = cos(mu*sqrt(tk))                      # theoretical result


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


#  untransformed plain Monte Carlo:

x = rnorm(N)                                 # using the same random numbers for all tk
Z = cosh( outer(x,sqrt(tk)) )
G = exp(mu^2/2) * cos(mu*x) 
G = outer(G,rep(1,nt))
EG0 = colSums(G*Z)/colSums(Z)                # untransformed result


# 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:

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