0. Introduction

This course project contains two topics of financial engineering models- Black-Scholes and Vasicek interest rate Model. Black-Scholes Model basically deals with option pricing methods, and Vasicek Model puts emphasis on interest rate term structure.

1. The Black-Scholes model for stock prices

Black-Scholes Model can apply to either call option or put option, and there are also some different versions of the model with few adjustments according to dividends (Black-Scholes-Merton Model) or being applied to Foreign Currency Exchange Option (Garman-Kohlhagen Model),….,etc. The following analysis only emphasizes on the computation of call option price without dividends.

i. Model briefing

The call option price in the Black-Scholes Model states as following:

\(P_{t0}\)=\(S_{0}\)N(\(\frac{log(S_{0}/c)+(\rho+\sigma^{2}/2)t_{0}}{\sigma\sqrt{t_{0}}}\)) - c\(e^{-\rho/t_{0}}\)N(\(\frac{log(S_{0}/c)+(\rho-\sigma^{2}/2)t_{0}}{\sigma\sqrt{t_{0}}})\)

,where N(.): cumulative normal distribution function, \(S_{0}\):the stock price at t=0, c: the strike price of the option, \(\rho\): interest rate, \(\sigma\): volatility

ii. Relative characteristics of the model (illustrated with plots)

The given parameter value is given as the following table.

##         Parameter Value
## 1     Stock price  1.00
## 2        Variance  0.02
## 3   Interest rate  0.03
## 4 Strike price(c)  1.00
set.seed(333)
StockPrice <- function(rho,t,sigma){
  mu= rho-sigma^2/2
  St= 1*exp(mu*t+sigma*rnorm(1,0,t))
  return(St)}
#Call price calculation function

call_price <- function(rho,sigma,c){
  #when t=0
  c0= max((1-c)*exp(-rho*0),0)
  C_price= c()
  for(i in 1:10){
    t=i
    S= StockPrice(rho,t,sigma)
    d1= (log(1/c)+(rho+sigma^2/2)*t)/(sigma*sqrt(t))
    d2= d1-sigma*sqrt(t)
    price= 1*pnorm(d1)-c*exp(-rho*t)*pnorm(d2)
    C_price[i]=price}
  C_price = append(C_price,c0, after=0)
  return(C_price)
}
#run with rho= 0.03, sigma= sqrt(0.02),c=1
c1 <- call_price(0.03,sqrt(0.02),1)
#run with rho= 0.02, sigma= sqrt(0.02), c=1
c2 <- call_price(0.02,sqrt(0.02),1)
#run with rho= 0.04, sigma= sqrt(0.02), c=1
c3 <- call_price(0.04,sqrt(0.02),1)


library(ggplot2)
g <- ggplot()
g+ geom_line(aes(x=0:10,y= c2, color= "black"))+
geom_line(aes(x=0:10, y=c1, color="blue"))+ geom_line(aes(x=0:10,y=c3,color= "red"))+ labs(title= "Call Price as time evolves(change with rho)", col= "rho")+ scale_color_manual(labels=c(0.02,0.03,0.04),values=c(black= "black",blue= "blue",red= "red"))+ xlab("Time")+ ylab("Call price")

#run with rho=0.03, sigma= sqrt(0.02),c=1
c1 <- call_price(0.03,sqrt(0.02),1)
#run with rho= 0.03, sigma= sqrt(0.01), c=1
c2 <- call_price(0.03,sqrt(0.01),1)
#run with rho= 0.03, sigma= sqrt(0.03), c=1
c3 <- call_price(0.03,sqrt(0.03),1)

library(ggplot2)
g <- ggplot()
g+ geom_line(aes(x=0:10,y= c2, color= "black"))+
geom_line(aes(x=0:10, y=c1, color="blue"))+ geom_line(aes(x=0:10,y=c3,color= "red"))+ labs(title= "Call Price as time evolves(change with volatility)", col= "sigma")+ scale_color_manual(labels=c("sqrt(0.01)","sqrt(0.02)","sqrt(0.03)"),values=c(black= "black",blue= "blue",red= "red"))+ xlab("Time")+ ylab("Call price")

The blue lines above are plotted according to \(S_{0}\) = 1, \(\sigma^2\)= 0.02, ρ = 0.03 and c = 1.

Since current stock price \(S_{0}\)=1, and the strike price c=1, the payoff of the call option on t=0 should be 0, thus the call price is 0 at t=0. From the formula (1i.) provided above, it is obvious that P_t tends to increase as t increases. Assume other things being equal, if the interest rate parameter \(\rho\) is shifted by a positive amount, the curve should move upward (the blue line shifts to the red line on the left plot above). When it comes to shifting the volatility, the higher the volatility, the higher the call price is.

#run with rho=0.03, sigma= sqrt(0.02),c=1
c1 <- call_price(0.03,sqrt(0.02),1)
#run with rho= 0.03, sigma= sqrt(0.02), c=0.5
c2 <- call_price(0.03,sqrt(0.02),0.5)
#run with rho= 0.03, sigma= sqrt(0.02), c=1.5
c3 <- call_price(0.03,sqrt(0.02),1.5)

library(ggplot2)
g <- ggplot()
g+ geom_line(aes(x=0:10,y= c2, color= "black"))+
geom_line(aes(x=0:10, y=c1, color="blue"))+ geom_line(aes(x=0:10,y=c3,color= "red"))+ labs(title= "Call Price as time evolves(change with strike price)", col= "strike price(c)")+ scale_color_manual(labels=c(0.5,1,1.5),values=c(black= "black",blue= "blue",red= "red"))+ xlab("Time")+ ylab("Call price")

Since when the strike price c=0.5, the holder of the option holder can attain profit=\(S_{0}\)-c=1-0.5=0.5 at time=0, the call price curve starts from 0.5 and increase as time evolves. If the strike price c= 1 or 1.5, the option holder will attain 0, so the call option price curve starts at 0 and increase as time evolves.

iii. Conclusion- Benefits and shortcomings

Black-Scholes Model provides a quick method to compute option price. Since it acts as a close solution, it is easy to calculate different prices by replacing parameters based on different market conditions. However, the model assumes the interest rate is always constant prior to the maturity, term structure should be taken into consideration if we want to compute the price of the option precisely.

2. The Vasicek model for interest rates

Vasicek model is used for constructing term structure of interest rate. In practical, interest rate won’t either always be flat or always moves in parellal shift. Term structure models are used to decide different interest rates according to different time.

i. Model briefing

Vasicek model defines the price, \(Q_{t}\), at time 0 of a bond paing one unit at time t as

\(Q_{t}\)= \(e^{-\left(\int_{0}^{t} R_{s}\; ds\right)}\)

\(R_{s}\) = \(e^{-\theta s}\)\(R_{0}\)+(1-\(e^{-\theta s}\))\(\mu\)+\(X_{s}\)

where {\(R_{s}\), s>0}: Ornstein-Uhlenbeck process, \(R_{0}\): initial spot rate, \(\mu\): long term mean level, \(\sigma\): speed of reversion(>0), \(X_{s}\): Ornstein-Uhlenbeck process with volatility>0 and reversion parameter \(\theta\)>0

ii. Relative characteristics of the model (illustrated with plots)

The following graphs are all plotted based on \(R_{0}\)=0.1,\(\ \theta\)=0.5, \(\mu\)=0.05, \(\sigma\)=0.02.

#OU process

rOU <- function(n,N,Delta,theta,sigma){
  times= (0:n)*Delta
  X= matrix(0, nrow=N, ncol= n+1)
  X[,1]= 0
  for(i in 1:n){
    x= X[,i]
    m= x*exp(-theta*Delta)
    v= sigma^2*(1-exp(-2*theta*Delta))/(2*theta)
    X[,i+1]= rnorm(N,m,sqrt(v))
  }
  return(X)
}

# vm
transf <- function(X,theta,mu,R0,s){exp(-theta*s)*R0+(1-exp(-theta*s))*mu+X}

result= c()
Rs <- function(n,N,Delta,theta,sigma,mu,R0,s){
  X= rOU(n,N,Delta,theta,sigma)
  j=1
  for(i in seq(0,s,Delta)){
    result[j]=transf(X[1,j],theta,mu,R0,i)
    j=j+1
  }
  return(result)
}
N=1
n=1000
Delta=0.01
R0=0.1
theta=0.5
mu=0.05
sigma= 0.02
s=10
times=(0:n)*Delta
y_vm = Rs(n,N,Delta,theta,sigma,mu,R0,s)

g <- ggplot()
g+ geom_line(aes(x=times, y=as.numeric(y_vm)))+labs(title="Realisations of Rs as time evolves- Graph(a.)")+ylab("Rs")

#Expected Rs 
E_Rs <- function(s,theta,R0,mu){   
  return(exp(-theta*s)*R0+(1-exp(-theta*s))*mu)}


plotY= list()
for(i in seq(0,s,Delta)){
  plotY= append(plotY,E_Rs(i,theta,R0,mu))
}

#theta=0.4
plotY1= list()
for(i in seq(0,s,Delta)){
  plotY1= append(plotY1,E_Rs(i,theta-0.1,R0,mu))
}

#theta= 0.6
plotY2= list()
for(i in seq(0,s,Delta)){
  plotY2= append(plotY2,E_Rs(i,theta+0.1,R0,mu))
}

g <- ggplot()
g+geom_line(aes(x=times, y=as.numeric(plotY1), color="black"))+geom_line(aes(x=times,y=as.numeric(plotY), color= "blue"))+ geom_line(aes(x=times,y=as.numeric(plotY2),color= "red"))+ labs(title= "E(Rs) change with theta- Graph(b.)",col= "theta")+ scale_color_manual(labels=c(0.4,0.5,0.6),values=c("blue","black","red"))+ xlab("Time")+ ylab("E(Rs)")

#change with mu
plotY= list()
for(i in seq(0,s,Delta)){
  plotY= append(plotY,E_Rs(i,theta,R0,mu))
}
#mu=0.04
plotY1= list()
for(i in seq(0,s,Delta)){
  plotY1= append(plotY1,E_Rs(i,theta,R0,mu-0.01))
}
#mu= 0.06
plotY2= list()
for(i in seq(0,s,Delta)){
  plotY2= append(plotY2,E_Rs(i,theta,R0,mu+0.01))
}

g <- ggplot()
g+geom_line(aes(x=times,y=as.numeric(plotY), color= "blue"))+geom_line(aes(x=times, y=as.numeric(plotY1), color="black"))+ geom_line(aes(x=times,y=as.numeric(plotY2),color= "red"))+ labs(title= "E(Rs) change with mu- Graph(c.)", col="mu")+ scale_color_manual(labels=c(0.04,0.05,0.06),values=c("blue", "black","red"))+ xlab("Time")+ ylab("E(Rs)")

# variance Rs
var_Rs <-  function(sigma,theta,Delta){
  return(sigma^2*(1-exp(-2*theta*Delta))/(2*theta))}

#change with sigma
plot_var= c()
for(i in seq(0,s,Delta)){
  plot_var = append(plot_var,var_Rs(sigma,theta,i))}

times=(0:n)*Delta

#sigma = 0.01
plot_var1= c()
for(i in seq(0,s,Delta)){
  plot_var1 = append(plot_var1,var_Rs(sigma-0.01,theta,i))}

#sigma =0.03
plot_var2= c()
for(i in seq(0,s,Delta)){
  plot_var2 = append(plot_var2,var_Rs(sigma+0.01,theta,i))}

times=(0:n)*Delta
g <- ggplot()
g+geom_line(aes(x=times, y=as.numeric(plot_var),col="blue"))+geom_line(aes(x=times, y=as.numeric(plot_var1),col="black"))+geom_line(aes(x=times, y=as.numeric(plot_var2),col="red"))+labs(title="Variance of Rs as time evolves(change with sigma)- Graph(d.)",col="sigma")+ylab("Variance of Rs")+scale_color_manual(labels=c(0.01,0.02,0.03),values=c("blue","black","red"))+ xlab("Time")+ ylab("E(Rs)")

#change with sigma


BP <- function(theta,mu,R0,z){
  bp <- function(t){
    exp(-theta*t)*R0+(1-exp(-theta*t))*mu+z}
  return(bp)}
  
z=0
count=0
Qt <- function(n,N,Delta,theta,sigma,mu,R0,s,z){
  y_bp=c()
  for(j in seq(0,s,Delta)){
    count= count+1
    z= rOU(n,N,Delta,theta,sigma)[count]
    intgral = integrate(BP(theta,mu,R0,z),lower=0,upper=j)$value
    y_bp[count]= exp(-1*intgral)
  }
  return(y_bp)
}


g <- ggplot()
g+ geom_line(aes(x=times,y=Qt(n,N,Delta,theta,sigma,mu,R0,s,z)))+ labs(title= "Bond Price related to time t- Graph(e.)")+ylab("Bond Price")

Qt2 <- function(n,N,Delta,theta,sigma,mu,R0,s,z){
  y_bp=c()
  x= rOU(n,N,Delta,theta,sigma)
  for(i in 1:N){
    z= x[i,ncol(rOU(n,N,Delta,theta,sigma))]
    intgral = integrate(BP(theta,mu,R0,z),lower=0,upper=s)$value
    y_bp[i]= exp(-1*intgral)}
  return(y_bp)
}
N=3000

g<- ggplot()
g+geom_histogram(aes(x=Qt2(n,N,Delta,theta,sigma,mu,R0,s,z)))+labs(title="Distribution of Qt- Graph(f.)")+xlab("Qt")+ylab("count")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Since Vasicek Model states that the interest rate \(R_{s}\) contains a stochastic factor \(X_{t}\), the realizations of \(R_{s}\) will also be a stochastic process starting from \(R_{0}\)=0.1 as Graph(a.) above.

As time increases, the expected value of \(R_{s}\) approaches \(R_{0}\) since E(\(X_{s}\))=0. Because the variance of \(R_{s}\) should turn out to be equal to var(\(X_{s}\))=$ (1-e^{-2s})$, it will approach to \(\frac{\sigma^2}{2\theta}\) as time increases. On Graph(b.), it starts from 0.1 and declines until 0.05 as time evolves. The variance starts from 0 and approaches to \(\frac{\sigma^2}{2\theta}=\frac{{0.02}^2}{2\times0.5}\)=0.0004, just as the black line (original parameter) shows on the Graph(d.). If the values of theta, mu and sigma are changed, the influences on E(\(R_{s}\)) and var(\(R_{s}\)) are illustrated as Graph(b.), Graph(c.) and Graph(d.) respectively.

The present value of a unit dollar paid on time t, denoted as Q_t ,is plotted on Graph(e.). Since the current one dollar does not have to be discounted, the graph starts from 1 and the value is discounted more and more as time increases. Due to Q_t contains  Ornstein-Uhlenbeck Process, {X_t}, the plot is consequently stochastic as it shows on Graph(e.). If t is fixed on 10, then the distribution of Q_t will look like Normal distribution (with 3000 times simulations), as it shows on Graph(f.).

iii. Conclusion- Benefits and shortcomings

Vasicek Model successfully captures the interest rate term structure. The benefit is that it takes the concept of mean reversion into the model, E(\(R_{s}\))=\(R_{0}\) as \(\ s\rightarrow\infty\), which is consistent to the observations and theories of mean-reverting in finance. The shortcoming of Vasicek Model is that the interest rate \(R_{s}\) can turn out to be 0 or negative, which violates the major theories in finance and the pricing conventions in fixed income products.