Assignment-2.utf8.md

Exercise 3

library(sde); library(dplyr)

## Loading required package: MASS

## Loading required package: stats4

## Loading required package: fda

## Warning: package 'fda' was built under R version 3.6.2

## Loading required package: splines

## Loading required package: Matrix

## Loading required package: fds

## Loading required package: rainbow

## Loading required package: pcaPP

## Warning: package 'pcaPP' was built under R version 3.6.2

## Loading required package: RCurl

## Warning: package 'RCurl' was built under R version 3.6.2

## Loading required package: deSolve

## Warning: package 'deSolve' was built under R version 3.6.2

## 
## Attaching package: 'fda'

## The following object is masked from 'package:graphics':
## 
##     matplot

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

## sde 2.0.15

## Companion package to the book

## 'Simulation and Inference for Stochastic Differential Equations With R Examples'

## Iacus, Springer NY, (2008)

## To check the errata corrige of the book, type vignette("sde.errata")

## Warning: package 'dplyr' was built under R version 3.6.2

## 
## Attaching package: 'dplyr'

## The following object is masked from 'package:MASS':
## 
##     select

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

#a)
#initialize varaibles 

mu = 0.15
sigma = 0.3
S0 = 100
T = 1 
nt = 3 
n = 250 
dt = T/n 

#generate data

epsilon = runif(n+1, 0, 1) #generate random deviates so the normal dist increments are different from each other 
epsilon = epsilon[order(epsilon)] #order dates in ascending order
X = matrix(rep(0,length(epsilon)*nt),nrow = nt)
X_star = matrix(rep(0,length(epsilon)*nt),nrow = nt) #X_star for when sigma = 0
t = seq(0, T, by = dt) 

delta_S = mu*S0*dt+sigma*epsilon*sqrt(dt)*S0
delta_S_perc = mu*dt+sigma*epsilon*sqrt(dt)

#loop

for(i in 1:nt) {
  X[i,] = sde::GBM(x = S0, r = mu, sigma = sigma, T = T, N = n)
}
for(i in 1:nt) {
  X_star[i,] = sde::GBM(x = S0, r = mu, sigma = 0, T = T, N = n)
}

ymax = max(X); ymin = min(X) #bounds for the prices

plot(epsilon,X[1,],type = "l", ylim = c(ymin,ymax), col = 1, xlab = "Time", ylab = "Price", main = "Brownian motion", cex.main = 1)

for(i in 2:nt){
  lines(epsilon,X[i,], type = 'l', ylim = c(ymin, ymax), col = i)
}
lines(t,X_star[1,], type = 'l', lty = 3, ylim = c(ymin, ymax), col = "blue") #line showing the trend corresponding to sigma = 0

X = as.data.frame(X)
X = t(X)
X_diff = diff(X) #The change in the price in the period

#b)
#sigma = 0.5%. Using the GBM function from the sde package
par(mfrow = c(3,1))
sigma2 = 0.005
X2 = matrix(rep(0,length(t)*nt),nrow = nt)
for(i in 1:nt) {
  X2[i,] = sde::GBM(x = S0, r = mu, sigma = sigma2, T = T, N = n)
}
plot(t,X2[1,],type = "l", col = 1, xlab = "Time", ylab = "Price", main = "Brownian motion with, sigma = .005", cex.main = 1)

#sigma = 20%. Using the GBM function from the sde package
sigma3 = 0.2
X3 = matrix(rep(0,length(t)*nt),nrow = nt)
for(i in 1:nt) {
  X3[i,] = sde::GBM(x = S0, r = mu, sigma = sigma3, T = T, N = n)
}
plot(t,X3[1,],type = "l", col = 1, xlab = "Time", ylab = "Price", main = "Brownian motion with sigma = .2", cex.main = 1)


#sigma = 50%. Using the GBM function from the sde package
sigma4 = 0.5
X4 = matrix(rep(0,length(t)*nt),nrow = nt)
for(i in 1:nt) {
  X4[i,] = sde::GBM(x = S0, r = mu, sigma = sigma4, T = T, N = n)
}
plot(t,X4[1,],type = "l", col = 1, xlab = "Time", ylab = "Price", main = "Brownian motion with sigma = .5", cex.main = 1)

Exercise 4

#a)
#Equation 3.4 (munk) states that the paths can be simulated by x_t_i = x_t_i-1 + μ(t_i-t_i-1)+ε_i*σ*sqrt(t_i-t_i-1)

rm(list = ls())
library(Sim.DiffProc)

## Warning: package 'Sim.DiffProc' was built under R version 3.6.2

## Package 'Sim.DiffProc', version 4.8
## browseVignettes('Sim.DiffProc') for more informations.

## 
## Attaching package: 'Sim.DiffProc'

## The following objects are masked from 'package:sde':
## 
##     BM, GBM

#initialize variables

mu = 10
sigma = 20
X0 = 100
T = 1 
nt = 3 
n = 250
dt = T/n

#generate data

epsilon = runif(n+1, 0, 1)
epsilon = epsilon[order(epsilon)]
X = matrix(rep(0,length(epsilon)*nt),nrow = nt)
X_star = matrix(rep(0,length(epsilon)*nt),nrow = nt)
t = seq(0, T, by = dt)

#loop

for(i in 1:nt) {
  X[i,] = ABM(x = X0, r = mu, sigma = sigma, T = T, N = n)
}
for(i in 1:nt) {
  X_star[i,] = ABM(x = X0, r = mu, sigma = 0, T = T, N = n)
}

ymax = max(X); ymin = min(X)

plot(epsilon,X[1,],type = "l", ylim = c(ymin,ymax), col = 1, xlab = "Time", ylab = "Price", main = "Brownian motion", cex.main = 1)

for(i in 2:nt){
  lines(epsilon,X[i,], type = 'l', ylim = c(ymin, ymax), col = i)
}
lines(t,X_star[1,], type = 'l', lty = 3, ylim = c(ymin, ymax), col = "blue")

#b) 

#scenerio 1, using the ABM function from Sim.diffProc package
par(mfrow = c(3,1))
sigma2 = 10
X2 = matrix(rep(0,length(t)*nt),nrow = nt)
for(i in 1:nt) {
  X2[i,] = ABM(x = X0, r = mu, sigma = sigma2, T = T, N = n)
}
plot(t,X2[1,],type = "l", col = 1, xlab = "Time", ylab = "Price", main = "Brownian motion with sigma = 10", cex.main = 1)



#scenerio 2
sigma3 = 20
X3 = matrix(rep(0,length(t)*nt),nrow = nt)
for(i in 1:nt) {
  X3[i,] = ABM(x = X0, r = mu, sigma = sigma3, T = T, N = n)
}
plot(t,X3[1,],type = "l", col = 1, xlab = "Time", ylab = "Price", main = "Brownian motion with sigma = 20", cex.main = 1)

#scenerio 3

sigma4 = 40
X4 = matrix(rep(0,length(t)*nt),nrow = nt)
for(i in 1:nt) {
  X4[i,] = ABM(x = X0, r = mu, sigma = sigma4, T = T, N = n)
}
plot(t,X4[1,],type = "l", col = 1, xlab = "Time", ylab = "Price", main = "Brownian sigma = 40", cex.main = 1)

Exercise 5

#a)
rm(list = ls())

#initialize variables
r0 = 0.05
mean = 0.05
k = 0.1
sigma = 0.05
 
nt = 3 
T = 1  
m = 250  
dt = T/m  

#generate data

t = seq(0, T, dt)
r = matrix(0,m+1,nt)
r[1,] = r0
 
#loop

for(j in 1:nt){
  for(i in 2:(m+1)){
    dr = k*(mean-r[i-1,j])*dt + sigma*sqrt(dt)*rnorm(1,0,1) #from the equation found on page 64 (munk) 
    r[i,j] = r[i-1,j] + dr
  }
} 
 
matplot(t, r[,1:nt], type = "l", lty = 1, main = "Ornstein-Uhlenbeck Process", ylab = "r_t", cex.main = 1)
abline(h = mean, col = "blue", lty = 2)
abline(h = 0, col = "red", lty = 2)

#b)
#scenerio 1
par(mfrow = c(3,1))
k2 = 0.1
for(j in 1:nt){
  for(i in 2:(m+1)){
    dr = k2*(mean-r[i-1,j])*dt + sigma*sqrt(dt)*sample(t, 1) #using sample() (as appose to rnorm()) to get the NID increments for the paths.
    r[i,j] = r[i-1,j] + dr
  }
}
matplot(t, r[,1], type = "l", lty = 1, main = "Interest Rate Paths with k=0.1", ylab = "r_t") 
abline(h = mean, col = "blue", lty = 2)

#scenerio 2
k3 = 10
for(j in 1:nt){
  for(i in 2:(m+1)){
    dr = k3*(mean-r[i-1,j])*dt + sigma*sqrt(dt)*sample(t, 1)
    r[i,j] = r[i-1,j] + dr
  }
}
matplot(t, r[,1], type = "l", lty = 1, main = "Interest Rate Paths with K=10", ylab = "r_t") 
abline(h = mean, col = "blue", lty = 2)

#scenerio 3
k4 = 200
for(j in 1:nt){
  for(i in 2:(m+1)){
    dr = k4*(mean-r[i-1,j])*dt + sigma*sqrt(dt)*sample(t, 1)
    r[i,j] = r[i-1,j] + dr
  }
}
matplot(t, r[,1], type = "l", lty = 1, main = "Interest Rate Paths with K=200", ylab = "r_t") 
abline(h = mean, col = "blue", lty = 2)

Exercise 6

#a)
rm(list = ls())

#initialize variables
r0 = 0.05
mean = 0.05
k = 0.1
sigma = 0.05
nt = 3 
T = 1   
m = 250   
dt = T/m  

#generate data

t = seq(0, T, dt) 
r = matrix(0,m+1,nt)
r[1,] = r0
 
#loop

for(j in 1:nt){
  for(i in 2:(m+1)){
    dr = k*(mean-r[i-1,j])*dt + sigma*sqrt(r[i-1,j])*sqrt(dt)*rnorm(1,0,1) #Equation from page 66 (Munk).
    r[i,j] = r[i-1,j] + dr
  }
} 
 
matplot(t, r[,1:nt], type = "l", lty = 1, main = "Square-root Process", ylab = "r_t") 
abline(h = mean, col = "blue", lty = 2)

#b)
#scenerio 1
par(mfrow = c(3,1))
k2 = 0.1
for(j in 1:nt){
  for(i in 2:(m+1)){
    dr = k2*(mean-r[i-1,j])*dt + sigma*sqrt(r[i-1,j])*sqrt(dt)*sample(t, 1)
    r[i,j] = r[i-1,j] + dr
  }
}
matplot(t, r[,1], type = "l", lty = 1, main = "Square-root Process with K=0.1", ylab = "r_t") 
abline(h = mean, col = "blue", lty = 2)

#scenerio 2
k3 = 10
for(j in 1:nt){
  for(i in 2:(m+1)){
    dr = k3*(mean-r[i-1,j])*dt + sigma*sqrt(r[i-1,j])*sqrt(dt)**sample(t, 1)
    r[i,j] = r[i-1,j] + dr
  }
}
matplot(t, r[,1], type = "l", lty = 1, main = "Square-root Process with K=10", ylab = "r_t") 
abline(h = mean, col = "blue", lty = 2)

#scenerio 3
k4 = 200
for(j in 1:nt){
  for(i in 2:(m+1)){
    dr = k4*(mean-r[i-1,j])*dt + sigma*sqrt(r[i-1,j])*sqrt(dt)**sample(t, 1)
    r[i,j] = r[i-1,j] + dr
  }
}
matplot(t, r[,1], type = "l", lty = 1, main = "Square-root Process with K=200", ylab = "r_t") 
abline(h = mean, col = "blue", lty = 2)