Stocks for the long run¶

To explore working with returns, we're going to simulate return data and look at how our portfolio would grow.

We will do analysis using relatively simple syntax common to many programming languages to make learning easier.

In [1]:

#r "nuget: FSharp.Stats, 0.5.0"
#r "nuget: Plotly.NET, 3.*"
#r "nuget: Plotly.NET.Interactive, 3.*"

open FSharp.Stats

open Plotly.NET

From 1/1871-1/2023, the US market return was annualized 7.5% with a 14% standard deviation (Robert Schiller data).

In [2]:

let rnorm = Distributions.Continuous.Normal.Init 0.075 0.14

for i = 0 to 5 do
    printfn $"{rnorm.Sample()}"

Out[2]:

0.11281501405599485-0.08408038022942653-0.13470733565870390.136710977869711760.106921953586262770.16073449907483076

Put sampled values in a list.

In [3]:

let returns = [ for i = 1 to 1000 do rnorm.Sample() ]

Plot the distribution of returns.

In [4]:

let hist = Chart.Histogram(returns)

In [ ]:

hist

Out[ ]:

In [6]:

let percentiles =
    [ for p in [0.0; 0.01; 0.05; 0.5; 0.95; 0.99; 1.0] do 
        let pctl = Quantile.compute p returns
        p, pctl ]

Out[6]:

val percentiles: (float * float) list =  [(0.0, -0.523598482); (0.01, -0.236118288); (0.05, -0.1491462588);   (0.5, 0.07794851357); (0.95, 0.3105842283); (0.99, 0.4150182609);   (1.0, 0.5562842221)]

In [7]:

for (p, pctl) in percentiles do
    printfn $"percentile %.2f{p} is {pctl}"

Out[7]:

percentile 0.00 is -0.5235984819519119percentile 0.01 is -0.23611828795794754percentile 0.05 is -0.14914625875278156percentile 0.50 is 0.07794851357062386percentile 0.95 is 0.31058422825831894percentile 0.99 is 0.41501826087627836percentile 1.00 is 0.5562842220620063

A draw of 10-years of returns.

In [8]:

let draw10y = [ for i = 1 to 10 do rnorm.Sample() ]

Out[8]:

val draw10y: float list =  [0.07090658771; 0.1290453357; 0.2513317343; 0.03518658973; 0.2223930737;   -0.06832573973; 0.1606886181; 0.07021140711; -0.03060338974;   -0.006223749925]

Compute the cumulative product of the returns.

In [9]:

let mutable cumprod = 1.0
for r in draw10y do
    cumprod <- cumprod * (1.0 + r)
    printfn $"r={round 3 r}"
    printfn $"    cumprod={round 3 cumprod}"

Out[9]:

r=0.071    cumprod=1.071r=0.129    cumprod=1.209r=0.251    cumprod=1.513r=0.035    cumprod=1.566r=0.222    cumprod=1.915r=-0.068    cumprod=1.784r=0.161    cumprod=2.07r=0.07    cumprod=2.216r=-0.031    cumprod=2.148r=-0.006    cumprod=2.135

Let's plot the cumulative return.

In [10]:

let draw10yCR = 
    let mutable year = 0
    let mutable cumprod = 1.0
    [ for r in draw10y do
        cumprod <- cumprod * (1.0 + r)
        year <- year + 1
        year, cumprod ]

Out[10]:

val draw10yCR: (int * float) list =  [(1, 1.070906588); (2, 1.209102088); (3, 1.512987812); (4, 1.566224694);   (5, 1.914542218); (6, 1.783729704); (7, 2.070354766); (8, 2.215717287);   (9, 2.147908827); (10, 2.13454078)]

Functional version of the cumulative return calculation.

In [11]:

((0, 1.0), draw10y)
||> List.scan (fun (year, cumprod) r -> 
    (year + 1, cumprod * (1.0 + r)))
|> List.tail

Out[11]:

val it: (int * float) list =  [(1, 1.070906588); (2, 1.209102088); (3, 1.512987812); (4, 1.566224694);   (5, 1.914542218); (6, 1.783729704); (7, 2.070354766); (8, 2.215717287);   (9, 2.147908827); (10, 2.13454078)]

Plot the cumulative return.

In [12]:

let draw10yCrPlot = Chart.Line(draw10yCR)

In [ ]:

draw10yCrPlot

Out[ ]:

Practice: Sample 20 years of returns and plot the cumulative return.

In [29]:

// Answer here

Let's simulate 1000 draws of 30 years of returns.

In [15]:

let draw1k = 
    [ for i = 1 to 1_000 do
        [ for y = 1 to 30 do rnorm.Sample() ] ]

let draw1kCR =
    [ for life in draw1k do
        let mutable accRet = 1.0
        let mutable year = 0
        [ for r in life do
            accRet <- accRet * (1.0 + r)
            year <- year+1
            year, accRet ] ]

Look at a few observtions from the first draw.

In [16]:

draw1kCR[0][0..5]

Out[16]:

val it: (int * float) list =  [(1, 1.176291981); (2, 1.270733795); (3, 1.499237526); (4, 1.923948337);   (5, 2.026791218); (6, 1.558203089)]

Plot the first draw

In [ ]:

Chart.Line(draw1kCR[0])

Out[ ]:

Plot the second draw

In [ ]:

Chart.Line(draw1kCR[1])

Out[ ]:

Plot the first and second together

In [ ]:

[ Chart.Line(draw1kCR[0])
  Chart.Line(draw1kCR[1]) ]
|> Chart.combine

Out[ ]:

Plot all 1000 draws together

In [ ]:

[ for x in draw1kCR do Chart.Line(x) ]
|> Chart.combine
|> Chart.withLayoutStyle(ShowLegend=false)

Out[ ]:

A more useful version may be to plot the values at the end of the last year.

In [21]:

let terminalValues = 
    [ for x in draw1kCR do 
        let (yr, ret) = x[x.Length-1]
        ret ]

In [ ]:

terminalValues |> Chart.Histogram

Out[ ]:

What is the chance we lose money?

In [23]:

let nLoseMoney =
    terminalValues
    |> List.filter (fun x -> x < 1.0)
    |> List.length
    |> float

let chanceLose = nLoseMoney / (float terminalValues.Length)    
printfn $"chance lose money=%.3f{chanceLose}"

Out[23]:

chance lose money=0.004

What is the chance we double our money?

In [24]:

let nDoubleMoney =
    terminalValues
    |> List.filter (fun x -> x >= 2.0)
    |> List.length
    |> float

let chanceDouble = nDoubleMoney / (float terminalValues.Length)    
printfn $"chance double money=%.3f{chanceDouble}"

Out[24]:

chance double money=0.945

Practice: What is the chance we earn at least a 10% compound rate of return per year?

In [29]:

// Answer here

Now let's say we have 1m EUR and we want to live off of 50k EUR per year. What is the chance that we can do that?

In [26]:

let expenses = 50_000.0
let initialWealth = 1_000_000.0
let wealthEvolution =
    [ for life in draw1k do
        let mutable wealth = initialWealth
        [ for r in life do
            // We'll take expenses out at the start of the year.
            if wealth > expenses then
                wealth <- (wealth - expenses) * (1.0 + r)
            else
                wealth <- 0.0
            wealth ] ]

let terminalWealth = [ for x in wealthEvolution do x[x.Length-1] ]

What's the chance that we don't have enough money?

In [27]:

let nBroke =
    terminalWealth
    |> List.filter (fun x -> x <= 0.0)
    |> List.length
    |> float

let chanceBroke = nBroke / (float terminalWealth.Length)    

printfn $"chance broke=%.3f{chanceBroke}"

Out[27]:

chance broke=0.158

Practice: What is our chance of going broke if the market's standard deviation is 20% per year?

In [29]:

// Answer here

Practice: What is our chance of going broke if the market's

return is 5% per year?

In [29]:

// Answer here