#--------------------------#
#   Loesungen UeBlatt 8    #
#--------------------------#


#  Aufgabe 1  #


# 1a)

f = function(x,y)
{
  w = sqrt(x^2+y^2)
  if(w==0)
  { res = 1 }                     # lim_{x->0} sin(x)/x = 1
  else
  { res = sin(w)/(w)  }
  return(res)
}

# 1b)

x = seq(-20,20,by=0.2)
y = x
nx = length(x)
ny = length(y)
z = matrix( 0 , nrow=nx , ncol=ny )

for(i in 1:nx)
{
  for(j in 1:ny)
  {
    z[i,j] = f(x[i],y[j])  
  }
}

persp(x,y,z)

require(plot3D)

persp3D(x,y,z)

# 1c) 

contour(x,y,z)

contour(x,y,z,nlevels=100)




#  Aufgabe 2  #


CompareBSPaths = function( mu, sigma, T, N, S0=100 )
{
  phi = rnorm(N,0,1)
  x = rep(0,N)                            # Brownian Motion
  S1 = rep(0,N)               
  S2 = rep(0,N)              
  dt = T/N
  tt = dt*(1:N)

  # Brownian Motion for all t_k:
  x = cumsum( sqrt(dt)*phi )           

  # BSPath mit Formel (3):
  S2 = S0*exp( (mu-sigma^2/2)*tt + sigma*x ) 
                                         
  # BSPath mit Formel (1):
  for( i in 1:N )
  {
    if(i==1)
      S1[i] =      S0 * ( 1 + mu*dt + sigma*sqrt(dt)*phi[i] )
    else
      S1[i] = S1[i-1] * ( 1 + mu*dt + sigma*sqrt(dt)*phi[i] )
  }

  result = list(S1,S2)
  names(result) = c("solution by iteration","closed form solution")

  info = "S1 in black, S2 in red"
  plot(tt, S1, ylim=c(0,300), main=info)
  lines(tt, S2, col="red")

  #return(result)   # we just look at the plot..
}

CompareBSPaths( mu=0.01, sigma=0.25 , T=5 , N=1000 )
CompareBSPaths( mu=0.01, sigma=0.25 , T=5 , N=100  )
CompareBSPaths( mu=0.01, sigma=0.25 , T=5 , N=10   )
CompareBSPaths( mu=0.01, sigma=0.25 , T=5 , N=10000)    # dt->0 <=> N->infty

CompareBSPaths( 0.15, 0.5 , 5 , 1000 )
CompareBSPaths( 0.15, 0.5 , 5 , 100 )
CompareBSPaths( 0.15, 0.5 , 5 , 10 )
CompareBSPaths( 0.15, 0.5 , 5 , 10000 )      # dt->0 <=> N->infty