In [1]:

```
%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$.

In [4]:

```
%%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

In [5]:

```
%%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.

In [8]:

```
%%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.

In [9]:

```
%%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 $$In [10]:

```
%%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.

In [11]:

```
%%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]$.

In [12]:

```
%%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

- $E(e_i^2) = Var(e_i) + E(e_i)^2 = 1 + 0 = 1$.
- $E(e_1 \cdot e_2) = Cov(e_1,e_2) + E(e_1)E(e_2) = 0 + 0 = 0$.
- $E(x_1) = E(e_1) = 0$.
- $E(x_2) = E(\rho \cdot e_1 + \sqrt{1-\rho^2}\cdot e_2)$.
- $Var(x_1) = E(x_1^2)-E(x_1)^2 = E(e_1^2) - 0 = 1-0 = 1$.
- $Var(x_2) = E(x_2^2)-E(x_2)^2 = E((\rho \cdot e_1 + \sqrt{1-\rho^2}\cdot e_2)^2)-0 = \rho^2 E(e_1^2)+(1-\rho^2)E(e_2^2) = 1$
- $Cov(x_1,x_2) = E(e_1 e_2) - E(e_1)E(e_2) = E[e_1(\rho \cdot e_1 + \sqrt{1-\rho^2}\cdot e_2)] = \rho E(e_1^2) + (1-\rho^2)E(e_1 \cdot e_2) = \rho$

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.

In [13]:

```
%%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)
```

We now compute the Cholesky decomposition of the covariance matrix.

In [14]:

```
%%R
L = t(chol(cv))
print(L)
```

The Cholesky decomposition is now used to generate multivariate random variables with the correct correlation.

In [15]:

```
%%R
e = matrix(rnorm(3*10000),10000,3)
x = t(L %*% t(e))
print(dim(x))
print(cov(x))
```

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.

In [16]:

```
%%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)
```

The expected return of the portfolio is

In [17]:

```
%%R
muP = t(w) %*% mu
print(muP)
```

[,1] [1,] 0.078

And the standard deviation of the portfolio is

In [18]:

```
%%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.

In [19]:

```
%%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.

In [20]:

```
%%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)
```

We can plot this to see the classic systematic risk plot. Systematic risk declines as the number of stocks in the portfolio increases.

In [21]:

```
%%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

- Increases when the relative return to it $(\mu - r_f {\bf 1})$ increases.
- Decreases when risk aversion increases.
- Decreases when riskiness of the assets increases as proxied for by $\Sigma$.

**Example**

In [22]:

```
%%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)
```

In [23]:

```
%%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 $$In [24]:

```
%%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.

In [25]:

```
%%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.

In [26]:

```
%%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

In [29]:

```
%%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)
}
```

In [30]:

```
%%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

In [ ]:

```
%%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)
}
```

In [32]:

```
%%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.

In [33]:

```
%%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.

In [34]:

```
%%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.

In [35]:

```
%%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)))
}
```

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.

In [36]:

```
%%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.

In [38]:

```
%%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.

In [39]:

```
nb_setup.images_hconcat(["DSTMAA_images/LSM1.png"], width=600)
```

Out[39]:

In [40]:

```
nb_setup.images_hconcat(["DSTMAA_images/LSM2.png"], width=600)
```

Out[40]: