%pylab inline
from ipypublish import nb_setup
%load_ext rpy2.ipython
Populating the interactive namespace from numpy and matplotlib
The law of motion for stocks is often based on a geometric Brownian motion, i.e.,
$$ dS(t) = \mu S(t) \; dt + \sigma S(t) \; dB(t), \quad S(0)=S_0. $$This is a stochastic differential equation (SDE), because it describes random movement of the stock $S(t)$. The coefficient $\mu$ determines the drift of the process, and $\sigma$ determines its volatility. Randomness is injected by Brownian motion $B(t)$. This is more general than a deterministic differential equation that is only a function of time, as with a bank account, whose accretion is based on the equation $dy(t) = r \;y(t) \;dt$, where $r$ is the risk-free rate of interest.
The solution to a SDE is not a deterministic function but a random function. In this case, the solution for time interval $h$ is known to be
$$ S(t+h) = S(t) \exp \left[\left(\mu-\frac{1}{2}\sigma^2 \right) h + \sigma B(h) \right] $$The presence of $B(h) \sim N(0,h)$ in the solution makes the function random. The presence of the exponential return makes the stock price lognormal. (Note: if random variable $x$ is normal, then $e^x$ is lognormal.)
Re-arranging, the stock return is
$$ R(t+h) = \ln\left( \frac{S(t+h)}{S(t)}\right) \sim N\left[ \left(\mu-\frac{1}{2}\sigma^2 \right) h, \sigma^2 h \right] $$Using properties of the lognormal distribution, the conditional mean of the stock price becomes
$$ E[S(t+h) | S(t)] = e^{\mu h} $$The following R code computes the annualized volatility $\sigma$.
%%R
library(quantmod)
getSymbols('MSFT',src='yahoo')
stkp = MSFT$MSFT.Close
rets = diff(log(stkp))[-1]
h = 1/252
sigma = sd(rets)/sqrt(h)
print(sigma)
[1] 0.2714769
The parameter $\mu$ is also easily estimated as
%%R
mu = mean(rets)/h+0.5*sigma^2
print(mu)
[1] 0.1597019
So the additional term $\frac{1}{2}\sigma^2$ does matter substantially.
It is easy to simulate a path of stock prices using a discrete form of the solution to the Geometric Brownian motion SDE.
$$ S(t+h) = S(t) \exp \left[\left(\mu-\frac{1}{2}\sigma^2 \right) h + \sigma \cdot e\; \sqrt{h} \right] $$Note that we replaced $B(h)$ with $e \sqrt{h}$, where $e \sim N(0,1)$. Both $B(h)$ and $e \sqrt{h}$ have mean zero and variance $h$. Knowing $S(t)$, we can simulate $S(t+h)$ by drawing $e$ from a standard normal distribution.
%%R
n = 252
s0 = 100
mu = 0.10
sig = 0.20
s = matrix(0,1,n+1)
h=1/n
s[1] = s0
for (j in 2:(n+1)) {
s[j]=s[j-1]*exp((mu-sig^2/2)*h
+sig*rnorm(1)*sqrt(h))
}
print(s[1:5])
print(s[(n-4):n])
plot(t(s),type="l",col="blue",xlab="Days",ylab="stock price"); grid(lwd=2)
[1] 100.0000 102.4160 103.3743 102.5718 102.4505 [1] 101.46903 99.49087 100.11088 100.49564 98.52290
The same logic may be used to generate multiple paths of stock prices, in a vectorized way as follows. In the following example we generate 3 paths. Because of the vectorization, the run time does not increase linearly with the number of paths, and in fact, hardly increases at all.
%%R
s = matrix(0,3,n+1)
s[,1] = s0
for (j in seq(2,(n+1))) {
s[,j]=s[,j-1]*exp((mu-sig^2/2)*h
+sig*matrix(rnorm(3),3,1)*sqrt(h))
}
ymin = min(s); ymax = max(s)
plot(t(s)[,1],ylim=c(ymin,ymax),type="l",xlab="Days",ylab="stock price"); grid(lwd=2)
lines(t(s)[,2],col="red",lty=2)
lines(t(s)[,3],col="blue",lty=3)
The plot code shows how to change the style of the path and its color.
If you generate many more paths, how can you find the probability of the stock ending up below a defined price? Can you do this directly from the discrete version of the Geometric Brownian motion process above?
To convert two independent random variables $(e_1,e_2) \sim N(0,1)$ into two correlated random variables $(x_1,x_2)$ with correlation $\rho$, use the following trannsformations.
$$ x_1 = e_1, \quad \quad x_2 = \rho \cdot x_1 + \sqrt{1-\rho^2} \cdot x_2 $$%%R
e = matrix(rnorm(20000),10000,2)
print(cor(e))
print(cor(e[,1],e[,2]))
[,1] [,2] [1,] 1.000000000 0.004768294 [2,] 0.004768294 1.000000000 [1] 0.004768294
We see that these are uncorrelated random variables. Now we create a pair of correlated variates using the formula above.
%%R
rho = 0.6
x1 = e[,1]
x2 = rho*e[,1]+sqrt(1-rho^2)*e[,2]
cor(x1,x2)
[1] 0.5963091
Algebraic and numerical check
Check algebraically that $E[x_i]=0, i=1,2$, $Var[x_i]=1, i=1,2$. Also check that $Cov[x_1,x_2]=\rho = Corr[x_1,x_2]$.
%%R
print(mean(x1))
print(mean(x2))
print(var(x1))
print(var(x2))
print(cov(x1,x2))
[1] -0.009192753 [1] 0.0001624843 [1] 0.9941886 [1] 1.019582 [1] 0.600367
We know that $E[e_i]=0, i=1,2$, $Var[e_i]=1, i=1,2$. Also
These are generated using Cholesky decomposition. We may write a covariance matrix in decomposed form, i.e., $\Sigma = L\; L'$, where $L$ is a lower triangular matrix. Alternatively we might have an upper triangular decomposition, where $U= L'$. Think of each component of the decomposition as a square-root of the covariance matrix.
%%R
#Original matrix
cv = matrix(c(0.01,0,0,0,0.04,0.02,0,0.02,0.16),3,3,byrow=TRUE)
print(cv)
#Let's enhance it
cv[1,2]=0.005
cv[2,1]=0.005
cv[1,3]=0.005
cv[3,1]=0.005
print(cv)
[,1] [,2] [,3]
[1,] 0.01 0.00 0.00
[2,] 0.00 0.04 0.02
[3,] 0.00 0.02 0.16
[,1] [,2] [,3]
[1,] 0.010 0.005 0.005
[2,] 0.005 0.040 0.020
[3,] 0.005 0.020 0.160
We now compute the Cholesky decomposition of the covariance matrix.
%%R
L = t(chol(cv))
print(L)
[,1] [,2] [,3] [1,] 0.10 0.00000000 0.0000000 [2,] 0.05 0.19364917 0.0000000 [3,] 0.05 0.09036961 0.3864367
The Cholesky decomposition is now used to generate multivariate random variables with the correct correlation.
%%R
e = matrix(rnorm(3*10000),10000,3)
x = t(L %*% t(e))
print(dim(x))
print(cov(x))
[1] 10000 3
[,1] [,2] [,3]
[1,] 0.009859271 0.004812116 0.005332888
[2,] 0.004812116 0.039911796 0.021392749
[3,] 0.005332888 0.021392749 0.162215532
In the last calculation, we confirmed that the simulated data has the same covariance matrix as the one that we generated correlated random variables from.
Let's enter a sample mean vector and covariance matrix and then using some sample weights, we will perform some basic matrix computations for portfolios to illustrate the use of R.
%%R
mu = matrix(c(0.01,0.05,0.15),3,1)
cv = matrix(c(0.01,0,0,0,0.04,0.02,
0,0.02,0.16),3,3)
print(mu)
print(cv)
w = matrix(c(0.3,0.3,0.4))
print(w)
[,1]
[1,] 0.01
[2,] 0.05
[3,] 0.15
[,1] [,2] [,3]
[1,] 0.01 0.00 0.00
[2,] 0.00 0.04 0.02
[3,] 0.00 0.02 0.16
[,1]
[1,] 0.3
[2,] 0.3
[3,] 0.4
The expected return of the portfolio is
%%R
muP = t(w) %*% mu
print(muP)
[,1] [1,] 0.078
And the standard deviation of the portfolio is
%%R
stdP = sqrt(t(w) %*% cv %*% w)
print(stdP)
[,1] [1,] 0.1868154
We thus generated the expected return and risk of the portfolio, i.e., the values $0.078$ and $0.187$, respectively.
We are interested in the risk of a portfolio, often measured by its variance. As we had seen in a previous chapter, as we increase $n$, the number of securities in the portfolio, the variance keeps dropping, and asymptotes to a level equal to the average covariance of all the assets. It is interesting to see what happens as $n$ increases through a very simple function in R that returns the standard deviation of the portfolio.
%%R
sigport = function(n,sig_i2,sig_ij) {
cv = matrix(sig_ij,n,n)
diag(cv) = sig_i2
w = matrix(1/n,n,1)
result = sqrt(t(w) %*% cv %*% w)
}
We now apply it to increasingly diversified portfolios.
%%R
n = seq(5,100,5)
risk_n = NULL
for (nn in n) {
risk_n = c(risk_n,sigport(nn,0.04,0.01))
}
print(risk_n)
[1] 0.1264911 0.1140175 0.1095445 0.1072381 0.1058301 0.1048809 0.1041976 [8] 0.1036822 0.1032796 0.1029563 0.1026911 0.1024695 0.1022817 0.1021204 [15] 0.1019804 0.1018577 0.1017494 0.1016530 0.1015667 0.1014889
We can plot this to see the classic systematic risk plot. Systematic risk declines as the number of stocks in the portfolio increases.
%%R
plot(n,risk_n,type="l",ylab="Portfolio Std Dev")
We will review the notation one more time. Assume that the risk free asset has return $r_f$. And we have $n$ risky assets, with mean returns $\mu_i, i=1...n$. We need to invest in optimal weights $w_i$ in each asset. Let $w = [w_1, \ldots, w_n]'$ be a column vector of portfolio weights. We define $\mu = [\mu_1, \ldots, \mu_n]'$ be the column vector of mean returns on each asset, and ${\bf 1} = [1,\ldots,1]'$ be a column vector of ones. Hence, the expected return on the portfolio will be
$$ E(R_p) = (1-w'{\bf 1}) r_f + w'\mu $$The variance of return on the portfolio will be
$$ Var(R_p) = w' \Sigma w $$where $\Sigma$ is the covariance matrix of returns on the portfolio. The objective function is a trade-off between return and risk, with $\beta$ modulating the balance between risk and return:
$$ U(R_p) = r_f + w'(\mu - r_f {\bf 1}) - \frac{\beta}{2} w'\Sigma w $$The f.o.c. becomes a system of equations now (not a single equation), since we differentiate by an entire vector $w$:
$$ \frac{dU}{dw'} = \mu - r_f {\bf 1} - \beta \Sigma w = {\bf 0} $$where the RHS is a vector of zeros of dimension $n$. Solving we have
$$ w = \frac{1}{\beta} \Sigma^{-1} (\mu - r_f {\bf 1}) $$Therefore, allocation to the risky assets
Example
%%R
n = 3
print(cv)
print(mu)
rf=0.005
beta = 4
wuns = matrix(1,n,1)
print(wuns)
w = (1/beta)*(solve(cv) %*% (mu-rf*wuns))
print(w)
w_in_rf = 1-sum(w)
print(w_in_rf)
[,1] [,2] [,3]
[1,] 0.01 0.00 0.00
[2,] 0.00 0.04 0.02
[3,] 0.00 0.02 0.16
[,1]
[1,] 0.01
[2,] 0.05
[3,] 0.15
[,1]
[1,] 1
[2,] 1
[3,] 1
[,1]
[1,] 0.1250000
[2,] 0.1791667
[3,] 0.2041667
[1] 0.4916667
%%R
beta = 3
w = (1/beta)*(solve(cv) %*% (mu-rf*wuns));
print(w)
beta = 2
w = (1/beta)*(solve(cv) %*% (mu-rf*wuns));
print(w)
[,1]
[1,] 0.1666667
[2,] 0.2388889
[3,] 0.2722222
[,1]
[1,] 0.2500000
[2,] 0.3583333
[3,] 0.4083333
Notice that the weights in stocks scales linearly with $\beta$. The relative proportions of the stocks themselves remains constant. Hence, $\beta$ modulates the proportions invested in a risk-free asset and a stock portfolio, in which stock proportions remain same. It is as if the stock versus bond decision can be taken separately from the decision about the composition of the stock portfolio. This is known as the two-fund separation property, i.e., first determine the proportions in the bond fund vs stock fund and the allocation within each fund can be handled subsequently.
Here are the equations for the Ornstein-Uhlenbeck (O-U) process, for mean-reverting interest rates, $r(t)$:
$$ dr(t) = k(\theta -r(t)) dt + \sigma dZ(t), \quad r(0)=r_0 $$where $k$ is the rate of mean reversion, $\theta$ is the long run mean of the process, and $\sigma$ is its volatility parameter. The solution to this stochastic differential equation is:
$$ r(t) = r(0) e^{-kt} + \theta (1-e^{-kt}) + \sigma \int_0^t e^{-k(t-s)}dZ(s) $$Note that the solution is a random function. We can use the solution to determine the conditional moments of the process as follows:
$$ E[r(t) | r(0)] = r(0) e^{-kt} + \theta (1-e^{-kt}) $$$$ Var[r(t)|r(0)] = \frac{\sigma^2}{2b}(1-e^{-2kt}) $$For simulation purposes, we discretize the solution as follows:
$$ r(t) = r(0) e^{-kt} + \theta (1-e^{-kt}) + \sigma e^{-0.5kt}\cdot z\sqrt{t} $$The O-U process permits negative interest rates. The CIR model handles this issue with a small modification:
$$ dr(t) = k(\theta -r(t)) dt + \sigma \sqrt{r} \cdot dZ(t), \quad r(0)=r_0 $$$$ r(t) = r(0) e^{-kt} + \theta (1-e^{-kt}) + \sigma e^{-0.5kt}\cdot z\sqrt{t\cdot r(t)} $$For the O-U process, a simple regression may be used to estimate its parameters. Discretizing the process leads to the following equation:
$$ r(t+h) - r(t) = k(\theta - r(t)) h + \sigma z \sqrt{h}, \quad z \sim N(0,1) $$We write $y(t)=r(t+h) - r(t)$, $\alpha = k \theta h$, $\beta = -kh$, and $e(t) = \sigma z(t) \sqrt{h}$. Then the equation above is written as
$$ y(t) = \alpha + \beta \cdot r(t) + e(t) $$which is a regression equation from which we can recover the coefficients of the model as follows:
\begin{eqnarray} k &=& -\frac{\beta}{h} \\ \theta &=& \frac{\alpha}{kh} = -\frac{\alpha}{\beta} \\ Var[e(t)] &=& \sigma^2 h \end{eqnarray}Therefore, also
$$ \sigma = \sqrt{\frac{1}{h}\cdot Var[e(t)]} = \frac{1}{\sqrt{h}}\cdot \sigma_e $$%%R
#FUNCTION TO READ IN CSV FILES FROM FRED
#Enter SeriesID as a text string
readFRED = function(SeriesID) {
url = paste("https://research.stlouisfed.org/fred2/series/",
SeriesID, "/downloaddata/",SeriesID,".csv",sep="")
data = readLines(url)
n = length(data)
data = data[2:n]
n = length(data)
df = matrix(0,n,2) #top line is header
for (j in 1:n) {
tmp = strsplit(data[j],",")
df[j,1] = tmp[[1]][1]
df[j,2] = tmp[[1]][2]
}
rate = as.numeric(df[,2])
idx = which(rate>0)
idx = setdiff(seq(1,n),idx)
rate[idx] = -99
date = df[,1]
df = data.frame(date,rate)
names(df)[2] = SeriesID
result = df
}
Download the one-year Treasury rates.
%%R
#Get the 1 yr rate
df = readFRED('DGS1')
idx = which(df$DGS1>0)
df = df[idx,]
print(dim(df))
plot(df, type='l')
grid()
[1] 14373 2
Run the estimation.
%%R
r = df[,2] #Change this to vary estimates
y = diff(r)
x = r[2:length(r)]
res = lm(y~x)
h = 1/260
a = res$coefficients[1]
b = res$coefficients[2]
sig_e = sd(res$residuals)
k = -b/h
theta = -a/b
sig = sig_e/sqrt(h)
print(c(k,theta,sig))
x (Intercept) -0.08004341 5.41126090 1.31564977
See: https://en.wikipedia.org/wiki/Black%E2%80%93Scholes_model
%%R
BS = function(S,K,T,v,rf,dv) {
d1 = (log(S/K) + (rf-dv+0.5*v^2)*T)/(v*sqrt(T))
d2 = d1 - v*sqrt(T)
bscall = S*exp(-dv*T)*pnorm(d1) - K*exp(-rf*T)*pnorm(d2)
bsput = -S*exp(-dv*T)*pnorm(-d1) + K*exp(-rf*T)*pnorm(-d2)
res = c(bscall,bsput)
}
%%R
S=50; K=51; T=1; v=0.20; rf=0.02; dv=0.0
cat("Call/Put price: ",BS(S,K,T,v,rf,dv))
Call/Put price: 3.987326 3.977459
%%R
simulateOption = function(m,n) {
h = T/n
stk = matrix(0,m,n+1)
e = matrix(rnorm(m*n),m,n)
stk[,1] = S
for (t in 2:(n+1)) {
stk[,t] = stk[,t-1]*exp((rf-dv-0.5*v^2)*h + v*sqrt(h)*e[,t-1])
}
c = mean(pmax(0,stk[,n+1]-K))*exp(-rf*T)
p = mean(pmax(0,K-stk[,n+1]))*exp(-rf*T)
res = c(c,p)
}
%%R
m = 50000
n = 50
c_vec = NULL
p_vec = NULL
for (j in 1:25) {
res = simulateOption(m,n)
c_vec = c(c_vec, res[1])
p_vec = c(p_vec, res[2])
}
cat("Calls: "); print(c(mean(c_vec), sd(c_vec)))
cat("Puts: "); print(c(mean(p_vec), sd(p_vec)))
Calls: [1] 3.98143171 0.04080088 Puts: [1] 3.98187113 0.02363918
The antithetic variable for a normal random variable $x$ is $x' = -x$. For a uniform random variable $u$, the antithetic is $u' = (1-u)$. In both cases, we are assured that the covariance between a variable and its antithetic is $Cov(x,x') \leq 0$. We can repeat the simulation model using the antithetic approach. The simulated estimates will have lower variance.
%%R
#Antithetic construction
antitheticOption = function(m,n) {
h = T/n
e = matrix(rnorm(m/2*n),m/2,n)
e = rbind(e,-e)
stk = matrix(0,m,n+1)
stk[,1] = S
for (t in 2:(n+1)) {
stk[,t] = stk[,t-1]*exp((rf-dv-0.5*v^2)*h + v*sqrt(h)*e[,t-1])
}
c = mean(pmax(0,stk[,n+1]-K))*exp(-rf*T)
p = mean(pmax(0,K-stk[,n+1]))*exp(-rf*T)
res = c(c,p)
}
c_vec = NULL
p_vec = NULL
for (j in 1:25) {
res = antitheticOption(m,n)
c_vec = c(c_vec, res[1])
p_vec = c(p_vec, res[2])
}
cat("Calls: "); print(c(mean(c_vec), sd(c_vec)))
cat("Puts: "); print(c(mean(p_vec), sd(p_vec)))
Calls: [1] 3.97960481 0.02303482 Puts: [1] 3.97210731 0.01450353
We see that the variance has come down with the antithetic approach.
This is an approach where a Monte Carlo pricing model $P$ is improved by benchmarking off a similar pricing model $Q$ for which a closed-form solution is known. The idea is to discover the simulation bias in model $Q$ and use it to correct the bias in model $P$.
Suppose the Monte Carlo estimate of the pricing problem is denoted $P$ and that of the analogous problem is denoted $Q$. The closed-form solution for the latter model is denoted $Q_0$. Error from simulation is proportional to $(Q-Q_0)$. We use this error to correct the price $P$ to a new price $P'$ as follows:
$$ P' = P + (Q_0 - Q) $$For this approach to work it must be that $Var(P') \leq Var(P)$, i.e., we achieve a reduction in pricing error. For this to be the case, the following condiiton must be satisfied:
$$ Var(P') = Var(P + (Q_0 - Q)) = Var(P-Q) = Var(P) + Var(Q) - 2 Cov(P,Q) $$This implies that for $Var(P')<Var(P)$ it must be that
$$ Cov(P,Q) > \frac{1}{2}Var(Q) $$In other words, model $Q$ must be closely related to model $P$ for this approach to be effective in reduction of pricing error.
A good example in option pricing for this approach is in the pricing of Asian options (i.e., options on the average stock price). There is no closed-form option pricing model for options on the arithmetic average price (model $P$) but there is for options on the geometric average price (model $Q$). The latter may be used to improve the pricing of the former.
American options are usually priced using backward recursion, whereas simulation runs forward in time. How might we adapt forward simulation to price American put options?
We note that an American put is exercised early at time $t$ when the stock price $S$ drops low enough that the gains from exercise $K-S(t)$ are greater than the expected future gains from exercise. The locus of points $B(t)$ where current gains from exercise are equal to future gains from exercise is known as the "early-exercise boundary" of the American put option. If we are given the boundary, then it is easy to simulate the stock price and every time it hits the boundary, we exercise and take profits. The average of the present value of these profits over many simulated paths is the price of the American put option.
But, we do not have the early-exercise boundary! So, in order to use simulation, we search over many boundaries till we find the one that gives the highest option value. We posit the following functional form for the boundary: $B(T-t) = K \exp(-\alpha (T-t))$. This simple functional form has a single parameter $\alpha$. By varying $\alpha$ we get a family of boundaries to search over. The program code is as follows.
%%R
simulateAmericanPut = function(m,n,alpha) {
h = T/n
stk = matrix(0,m,n+1)
stk[,1] = S
for (t in 2:(n+1)) {
stk[,t] = stk[,t-1]*exp((rf-dv-0.5*v^2)*h + v*sqrt(h)*e[,t-1])
}
total = 0
for (t in 2:(n+1)) {
idx = which(stk[,t] < K*exp(-alpha*(T-(t-1)*h)))
if (length(idx)>0) {
total = total + sum(K-stk[idx,t])*exp(-rf*(t-1)*h)
stk = stk[-idx,]
}
}
ap = (total/m)
}
Let's run this for different $\alpha$ to find the optimal early-exercise boundary and the price of the American put.
%%R
e = matrix(rnorm(m*n),m,n)
for (alpha in seq(0.10,1.00,0.05)) {
print(c(alpha, simulateAmericanPut(m,n,alpha)))
}
[1] 0.100000 3.200471 [1] 0.150000 3.575733 [1] 0.200000 3.785341 [1] 0.250000 3.896638 [1] 0.300000 3.958082 [1] 0.350000 3.997789 [1] 0.400000 4.021776 [1] 0.450000 4.033384 [1] 0.500000 4.039299 [1] 0.550000 4.039897 [1] 0.600000 4.039012 [1] 0.650000 4.034332 [1] 0.700000 4.033307 [1] 0.750000 4.033123 [1] 0.800000 4.035838 [1] 0.850000 4.028741 [1] 0.900000 4.026945 [1] 0.950000 4.028077 [1] 1.000000 4.029204
Tha value here is greater than that of the European put option.
We can reverify these results using a binomial tree approach to pricing the American put.
%%R
CRR = function(S0,K,T,v,rf,n) {
h = T/n
u = exp(v*sqrt(h))
d = 1/u
R = exp(rf*h)
q = (R-d)/(u-d)
#Stock grid
S = matrix(0,n+1,n+1)
P = S
S[1,1] = S0
for (t in 2:(n+1)) {
S[1,t] = S[1,t-1]*u
for (j in 2:t) {
S[j,t] = S[j-1,t-1]*d
}
}
#Terminal Value
P[,n+1] = pmax(0,K-S[,n+1])
#Backward recursion
for (t in seq(n,1,-1)) {
P[1:t,t] = (q*P[1:t,t+1]+(1-q)*P[2:(t+1),t+1])*exp(-rf*h)
idx = which(K-S[1:t,t] > P[1:t,t])
P[idx,t] = K-S[idx,t]
}
res = P[1,1]
}
Run the program to get the correct American put price.
%%R
print(CRR(S,K,T,v,rf,50))
[1] 4.100883
As we see, the price is close to the one obtained by simulation.
In an important paper, Francis Longstaff and Eduardo Schwartz developed a simple and interesting approach to pricing American options using regressions. We present the example from their paper here, see @doi:10.1093/rfs/14.1.113.
nb_setup.images_hconcat(["DSTMAA_images/LSM1.png"], width=600)
nb_setup.images_hconcat(["DSTMAA_images/LSM2.png"], width=600)
nb_setup.images_hconcat(["DSTMAA_images/LSM3.png"], width=600)
The entire paper can be programmed as follows.
%%R
#LSM MODEL
#FUNCTION TO SIMULATE A SET OF m PATHS OF STOCK PRICE FOR n PERIODS
stock_grid = function(S0,rf,v,h,n,m) {
s = matrix(0,m,n+1) #Blank stock matrix
e = matrix(rnorm(m*n),m,n) #Random normal set
s[,1] = S0
for (j in 1:n) {
s[,j+1] = s[,j]*exp((rf-0.5*v^2)*h + v*sqrt(h)*e[,j]) #Geometric Brownian motion
}
result = s
}
#LSM MODEL
LSM = function(S0,K,T,v,rf,h,n,m) {
# Generate stock grid
stock = stock_grid(S0,rf,v,h,n,m)
# Initialize option grid
option = matrix(0,m,n+1);
# Compute value of option at maturity
option[,n+1] = pmax(0, K - stock[,n+1])
# Backward recursion
for (i in seq(n, 2, -1)) {
# regression step
index = which(K - stock[, i] > 0)
# continuation value of option
W = exp(-rf * h) * option[, i + 1];
# the Y vector of continuation values for regression
Y = W[index];
S = stock[, i];
# the stock price (S) vector for regression
S = S[index];
# matrix for 3rd order polynomial regression
S = matrix(c(S, S^2, S^3), length(S), 3);
beta = lm(Y ~ S);
# continuation value step
earlyIndex = matrix(0, m, 1);
earlyIndex[index] = 1;
S = stock[, i];
S = matrix(c(S, S^2, S^3), length(S), 3);
predictedY = S %*% beta$coefficients[2:4] + beta$coefficients[1];
contVal = predictedY * earlyIndex;
earlyIndex2 = matrix(0, m, 1);
index = which(K - stock[, i] > contVal);
earlyIndex2[index] = 1;
earlyIndex = earlyIndex * earlyIndex2;
option[, i] = earlyIndex * (K - stock[, i]) + (1 - earlyIndex) * W;
}
# Complete option value at initial date
option[, 1] = option[, 2] * exp(-rf * h);
result = max(mean(option[, 1]), K - S0);
}
%%R
h = T/n
putval = LSM(S,K,T,v,rf,h,n,m)
print(c("Put value = ", putval))
[1] "Put value = " "4.07854140591467"
The general approach for Monte Carlo approaches for option pricing may be referenced here: https://en.wikipedia.org/wiki/Monte_Carlo_methods_for_option_pricing