In [1]:
using Dates, DelimitedFiles, LinearAlgebra, Roots, Distributions, StatsBase
include("printmat.jl")
include("printTable.jl")
Out[1]:
In [2]:
using Plots
#pyplot(size=(600,400)) #use pyplot or gr
gr(size=(480,320))
default(fmt = :svg)
In [3]:
x = readdlm("Data/MyData.csv",',',skipstart=1) #monthly return data
ym = round.(Int,x[:,1]) #yearmonth, like 200712
Rme = x[:,2] #market excess return
Rf = x[:,3] #interest rate
R = x[:,4] #return small growth stocks
Re = R - Rf #excess returns
T = size(Rme,1)
dN = Date.(string.(ym),"yyyymm") #convert to string and then Julia Date
printmat(dN[1:4])
In [4]:
x = [ones(T) Rme] #regressors
y = copy(Re) #to get standard OLS notation
b = x\y #OLS
u = y - x*b #residuals
covb = inv(x'x)*var(u) #cov(b), see any textbook
stdb = sqrt.(diag(covb)) #std(b)
R2a = 1 - var(u)/var(y)
printTable([b stdb b./stdb],["coeff","std","t-stat"],["α","β"])
printlnPs("R2: ",R2a)
printlnPs("no. of observations: ",T)
In [5]:
plags = 1:5
xCorr = autocor(Re,plags) #using the StatsBase package
println("Autocorrelations (different lags) of the excess returns in Re")
printTable([xCorr sqrt(T)*xCorr],["autocorr","t-stat"],string.(plags),cell00="lag")
The next cell implements a very simple momentum trading strategy.
If $R_{t-1}^{e}\ge0$, then we hold the market index and shorten the riskfree from $t-1$ to $t$. This means that we will earn $R_{t}^{e}$.
Instead, if $R_{t-1}^{e}<0$, then we do the opposite. This means that we will earn $-R_{t}^{e}$.
This simple strategy could be coded without using a loop, but "vectorization" does not speed up much.
In [6]:
(w,Rp) = (fill(NaN,T),fill(NaN,T))
for t = 2:T
w[t] = (Re[t-1] < 0)*(-1) + (Re[t-1] >= 0)*1 #w is -1 or 1
Rp[t] = w[t]*Re[t]
end
μ = [mean(Rp[2:end]) mean(Re[2:end])]
σ = [std(Rp[2:end]) std(Re[2:end])]
printlnPs("The annualized mean excess return of the strategy and a passive portfolio are: ",μ*12)
printlnPs("The annualized Sharpe ratios are: ",sqrt(12)*μ./σ)
The next cell constructs an simple estimate of $\sigma_t^2$ as a backward looking moving average (the RiskMetrics approach):
$\sigma_t^2 = \lambda \sigma_{t-1}^2 + (1-\lambda) (R_{t-1} -\mu)^2$, where $\mu$ is the average return (for all data).
Then, we calculate the 95% VaR by assuming a $N(\mu,\sigma_t^2)$ distribution:
$\textrm{VaR}_{t} = - (\mu-1.64\sigma_t)$.
If the model is correct, then $-R_t > \text{VaR}_{t}$ should only happen 5% of the times.
In [7]:
μ = mean(Rme)
λ = 0.95 #weight on old volatility
σ² = fill(var(Rme),T) #RiskMetrics approach to estimate variance
for t = 2:T
σ²[t] = λ*σ²[t-1] + (1-λ)*(Rme[t-1]-μ)^2
end
VaR95 = -(μ .- 1.64*sqrt.(σ²)) #VaR at 95% level
println(" ")
In [8]:
xTicksLoc = Dates.value.([Date(1980);Date(1990);Date(2000);Date(2010)])
xTicksLab = ["1980";"1990";"2000";"2010"] #crude way of getting tick marks right
plot( dN,VaR95,
color = :blue,
legend = false,
xticks = (xTicksLoc,xTicksLab),
ylim = (0,11),
title = "1-month Value at Risk (95%)",
ylabel = "%",
annotation = (Dates.value(Date(1982)),1,text("(for US equity market)",8,:left)) )
Out[8]:
Let $S$ be the the current spot price of an asset and $y$ be the interest rate.
The Black-Scholes formula for a European call option with strike price $K$ and time to expiration $m$ is
$C =S\Phi(d_{1}) -e^{-ym}K\Phi(d_{2})$, where
$d_{1} =\frac{\ln(S/K)+(y+\sigma^{2}/2)m}{\sigma\sqrt{m}} \ \text{ and } \ d_{2}=d_{1}-\sigma\sqrt{m}$
and where $\Phi(d)$ denotes the probability of $x\leq d$ when $x$ has an $N(0,1)$ distribution. All variables except the volatility ($\sigma$) are directly observable.
In [9]:
Φ(x) = cdf(Normal(0,1),x) #a one-line function, Pr(z<=x) for N(0,1)
"""
Calculate Black-Scholes european call option price
"""
function OptionBlackSPs(S,K,m,y,σ)
d1 = ( log(S/K) + (y+1/2*σ^2)*m ) / (σ*sqrt(m))
d2 = d1 - σ*sqrt(m)
c = S*Φ(d1) - exp(-y*m)*K*Φ(d2)
return c
end
Out[9]:
In [10]:
σ = 0.4
c1 = OptionBlackSPs(10,10,0.5,0.1,σ)
printlnPs("\n","call price according to Black-Scholes: ",c1)
K = range(7,stop=13,length=51)
c = OptionBlackSPs.(10,K,0.5,0.1,σ);
In [11]:
plot( K,c,
color = :red,
legend = false,
title = "Black-Scholes call option price",
xlabel = "strike price",
ylabel = "option price" )
Out[11]:
is the $\sigma$ value that makes the Black-Scholes equation give the same option price as observed on the market. It is often interpreted as the "market uncertainty."
The next cell uses the call option price calculated above as the market price. The implied volatility should then equal the volatility used above (this is a way to check your coding).
The next few cells instead use some data on options on German government bonds.
In [12]:
#solve for implied vol
iv = find_zero(σ->OptionBlackSPs(10,10,0.5,0.1,σ)-c1,(0.00001,5))
printlnPs("Implied volatility: ",iv,", compare with: $σ")
In [13]:
# LIFFE Bunds option data, trade date April 6, 1994
K = [ #strike prices; Mx1 vector
92.00; 94.00; 94.50; 95.00; 95.50; 96.00; 96.50; 97.00;
97.50; 98.00; 98.50; 99.00; 99.50; 100.0; 100.5; 101.0;
101.5; 102.0; 102.5; 103.0; 103.5 ];
C = [ #call prices; Mx1 vector
5.13; 3.25; 2.83; 2.40; 2.00; 1.64; 1.31; 1.02;
0.770; 0.570; 0.400; 0.280; 0.190; 0.130; 0.0800; 0.0500;
0.0400; 0.0300; 0.0200; 0.0100; 0.0100 ];
S = 97.05 #spot price
m = 48/365 #time to expiry in years
y = 0.0 #Interest rate: LIFFE=>no discounting
N = length(K)
Out[13]:
In [14]:
iv = fill(NaN,N) #looping over strikes
for i = 1:N
iv[i] = find_zero(sigma->OptionBlackSPs(S,K[i],m,y,sigma)-C[i],(0.00001,5))
end
println("Strike and iv for data: ")
printmat([K iv])
In [15]:
plot( K,iv,
color = :red,
legend =false,
title = "Implied volatility",
xlabel = "strike price",
annotation = (98,0.09,text("Bunds options April 6, 1994",8,:left)) )
Out[15]:
Given a vector of average returns ($\mu$) and a variance-covariance matrix ($\Sigma$), the mean-variance frontier shows the lowest possible portfolio uncertainty for a given expected portfolio return (denoted $\mu\text{star}$ below).
It is thus the solution to a quadratic minimization problem. The cells below will use the explicit (matrix) formulas for this solution, but we often have to resort to numerical methods when there are portfolio restrictions.
It is typically plotted with the portfolio standard deviation on the horizontal axis and the portfolio expected return on the vertical axis.
We calculate and plot two different mean-variance frontiers: (1) when we only consider risky assets; (2) when we also consider a risk-free asset.
In [16]:
μ = [11.5; 9.5; 6]/100 #expected returns
Σ = [166 34 58; #covariance matrix
34 64 4;
58 4 100]/100^2
Rf = 0.03 #riskfree return (an interest rate)
println("μ: ")
printmat(μ)
println("Σ: ")
printmat(Σ)
println("Rf: ")
printmat(Rf)
In [17]:
function MVCalc(μstar,μ,Σ) #calculate the std of a portfolio on MVF of risky assets
n = length(μ) #μstar is a scalar, μ a vector and Σ a matrix
oneV = ones(n)
Σ_1 = inv(Σ)
A = μ'Σ_1*μ
B = μ'Σ_1*oneV
C = oneV'Σ_1*oneV
λ = (C*μstar - B)/(A*C-B^2)
δ = (A-B*μstar)/(A*C-B^2)
w = Σ_1 *(μ*λ + δ*oneV)
StdRp = sqrt(w'Σ*w)
return StdRp,w
end
function MVCalcRf(μstar,μ,Σ,Rf) #calculates the std of a portfolio on MVF of (Risky,Riskfree)
n = length(μ)
μe = μ .- Rf #expected excess returns
Σ_1 = inv(Σ)
w = (μstar - Rf)/(μe'Σ_1*μe) * Σ_1 *μe
StdRp = sqrt(w'Σ*w)
return StdRp,w
end
Out[17]:
In [18]:
μstar = range(Rf,stop=0.15,length=201)
L = length(μstar)
StdRp = fill(NaN,L) #risky assets only
for i = 1:L #loop over different expected portfolio returns
StdRp[i] = MVCalc(μstar[i],μ,Σ)[1] #[1] to get the first output
end
In [19]:
plot( StdRp*100,μstar*100,
linecolor = :red,
xlim = (0,15),
ylim = (0,15),
label = "MVF",
legend = :topleft,
title = "MVF, only risky assets",
xlabel = "Std(Rp), %",
ylabel = "ERp, %" )
scatter!(sqrt.(diag(Σ))*100,μ*100,color=:red,label="assets")
Out[19]:
In [20]:
StdRpRf = fill(NaN,L) #with riskfree too
for i = 1:L
StdRpRf[i] = MVCalcRf(μstar[i],μ,Σ,Rf)[1]
end
In [21]:
plot( [StdRp StdRpRf]*100,[μstar μstar]*100,
legend = nothing,
linestyle = [:solid :dash],
linecolor = [:red :blue],
xlim = (0,15),
ylim = (0,15),
title = "MVF, risky and riskfree assets",
xlabel = "Std(Rp), %",
ylabel = "ERp, %" )
scatter!(sqrt.(diag(Σ))*100,μ*100,color=:red)
Out[21]:
The code below solves (numerically) the following minimization problem
$\min \text{Var}(R_p) \: \text{ s.t. } \: \text{E}R_p = \mu^*$,
and where we also require $w_i\ge 0$ and $\sum_{i=1}^{n}w_{i}=1$.
To solve this, we use the JuMP package (version >= 0.19) together with Ipopt.
In [22]:
using JuMP, Ipopt
In [23]:
function MeanVarNoSSPs(μ,Σ,μstar) #MV with no short-sales, numerical minimization
n = length(μ)
if minimum(μ) <= μstar <= maximum(μ) #try only if feasible
model = Model(with_optimizer(Ipopt.Optimizer, print_level = 1))
@variable(model,w[i=1:n] >= 0.0) #no short sales
@objective(model,Min,w'*Σ*w) #minimize portfolio variance
@constraint(model,sum(w) == 1.0) #w sums to 1
@constraint(model,w'*μ == μstar) #mean equals μstar
optimize!(model)
if has_values(model)
w_p = value.(w)
StdRp = sqrt(w_p'Σ*w_p)
end
else
(w_p,StdRp) = (NaN,NaN)
end
return StdRp,w_p
end
Out[23]:
In [24]:
Std_no_ss = fill(NaN,length(μstar))
for i = 1:length(μstar) #risky assets only, no short sales
Std_no_ss[i] = MeanVarNoSSPs(μ,Σ,μstar[i])[1]
end
In [25]:
plot( [StdRp Std_no_ss]*100,[μstar μstar]*100,
linecolor = [:red :green],
linestyle = [:solid :dash],
linewidth = 2,
label = ["no constraints" "no short sales"],
xlim = (0,15),
ylim = (0,15),
legend = :topleft,
title = "MVF (with/without constraints)",
xlabel = "Std(Rp), %",
ylabel = "ERp, %" )
scatter!(sqrt.(diag(Σ))*100,μ*100,color=:red,label="assets")
Out[25]:
In [ ]: