Introduction

This notebook is an implementation of Jon Danielsson's Financial Risk Forecasting (Wiley, 2011) in Julia v1.7.3, with annotations and introductory examples. The introductory examples (Appendix) are similar to Appendix B/C in the original book, with an emphasis on the differences between R/MATLAB and Julia.

Bullet point numbers correspond to the R/MATLAB Listing numbers in the original book, referred to henceforth as FRF.

More details can be found at the book website: https://www.financialriskforecasting.com/

Last updated: June 2022

Copyright 2011-2022 Jon Danielsson. This code is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. This code is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. The GNU General Public License is available at: https://www.gnu.org/licenses/.


Appendix: An Introduction to Julia

Created in Julia 1.7.3 (June 2022)

  • J.1: Entering and Printing Data
  • J.2: Vectors, Matrices and Sequences
  • J.3: Importing Data (to be updated)
  • J.4: Basic Summary Statistics
  • J.5: Calculating Moments
  • J.6: Basic Matrix Operations
  • J.7: Statistical Distributions
  • J.8: Statistical Tests
  • J.9: Time Series
  • J.10: Loops and Functions
  • J.11: Basic Graphs
  • J.12: Miscellaneous Useful Functions
In [1]:
# Entering and Printing Data in Julia
# Listing J.1
# Last updated June 2018
#
#

x = 10             # assign x the value 10
println(x)         # print x

## println() puts next output on new line, while print() doesn't
10
In [2]:
# Vectors, Matrices and Sequences in Julia
# Listing J.2
# Last updated July 2020
#
#

y = [1,3,5,7,9]     # lists in square brackets are stored as arrays

println(y)

println(y[3])       # calling 3rd element (Julia indices start at 1)

println(size(y))    # size of y

println(length(y))  # as expected, y has length 5

v = fill!(Matrix{Float64}(undef, 2,3),NaN) # 2x3 Float64 matrix of NaNs - computationally better
v = fill(NaN, (2,3))                       # 2x3 Float64 matrix of NaNs - direct

println(v)          # Julia prints matrices in a single line

println(size(v))    # as expected, v is size (2,3)

w = repeat([1,2,3]', outer = [3,2]) # repeats matrix thrice by rows, twice by columns

println(w)

s = 1:10            # s is an sequence which one can loop across

println(collect(s)) # return sequence elements as an array
[1, 3, 5, 7, 9]
5
(5,)
5
[NaN NaN NaN; NaN NaN NaN]
(2, 3)
[1 2 3 1 2 3; 1 2 3 1 2 3; 1 2 3 1 2 3]
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
In [3]:
# Basic Summary Statistics in Julia
# Listing J.3
# Last updated July 2020
#
#

y = [3.14,15,9.26,5]

using Statistics;                # load package needed
println("sum: ", sum(y))         # return sum of all elements of y
println("product: ", prod(y))    # return product of all elements of y
println("max: ", maximum(y))     # return maximum value of y
println("min: ", minimum(y))     # return minimum value of y
println("mean: ", mean(y))       # arithmetic mean
println("median: ", median(y))   # median
println("variance: ", var(y))    # variance
println("cov_matrix: ", cov(y))  # covar matrix = variance for single vector
println("cor_matrix: ", cor(y))  # corr matrix = [1] for single vector
println(sort(y))                 # sorts y in ascending order
println(log.(y))                 # natural log, note . denotes elementwise operation
sum: 32.4
product: 2180.73
max: 15.0
min: 3.14
mean: 8.1
median: 7.13
variance: 27.722400000000004
cov_matrix: 27.722400000000004
cor_matrix: 1.0
[3.14, 5.0, 9.26, 15.0]
[1.144222799920162, 2.70805020110221, 2.225704048658088, 1.6094379124341003]
In [4]:
# Calculating Moments in Julia
# Listing J.4
# Last updated June 2018
#
#

using StatsBase;

println("mean: ", mean(y))          # mean
println("variance: ", var(y))       # variance
println("std dev: ", std(y))        # unbiased standard deviation
println("skewness: ", skewness(y))  # skewness
println("kurtosis: ", kurtosis(y))  # EXCESS kurtosis (note the different default)
mean: 8.1
variance: 27.722400000000004
std dev: 5.265206548655049
skewness: 0.47004939528995604
kurtosis: -1.2846861061904815
In [5]:
# Basic Matrix Operations in Julia
# Listing J.5
# Last updated June 2018
#
#

z = Matrix([[1 2];[3 4]])   # z is a 2 x 2 matrix
x = Matrix([1 2])           # x is a 1 x 2 matrix

## Note: z * x is undefined since the two matrices are not conformable

println(z * x')             # this evaluates to a 2 x 1 matrix

b = vcat(z,x)               # "stacking" z and x vertically
c = hcat(z,x')              # "stacking" z and x' horizontally 

## Note: dimensions must match along the combining axis
[5; 11;;]
Out[5]:
2×3 Matrix{Int64}:
 1  2  1
 3  4  2
In [6]:
# Statistical Distributions in Julia
# Listing J.6
# Last updated July 2020
#
#

## Julia has a wide range of functions contained in the package Distributions.jl
## Vectorized versions of the functions are used here as they are relevant for FRF

using Distributions;

q = collect((-3:1:3))              # specify a set of values

p = collect((0.1:0.1:0.9))         # specify a set of probabilities

println(quantile.(Normal(0,1),p))  # element-wise inverse Normal quantile

println(cdf.(TDist(4), q))         # element-wise cdf calculation under Student-t(4)

println(pdf.(Chisq(2), q))         # element-wise pdf calculation under Chisq(2)

## Similar syntax for other dists, e.g. Bernoulli(p), Binomial(n,p), Poisson(λ)

## For full list of supported distributions, see Distributions.jl documentation

## One can also obtain pseudorandom samples from distributions using rand()

x = rand(TDist(5), 100)            # Sampling 100 times from TDist with 5 df

y = rand(Normal(0,1), 50)          # Sampling 50 times from a standard normal 

## Given data, we obtain MLE estimates of parameters with fit_mle():

fit_mle(Normal, x)                 # Fitting x to normal dist     

## Some distributions like the Student-t cannot be fitted yet (as of July 2020)
## Supported dists: https://juliastats.org/Distributions.jl/stable/fit/
[-1.2815515655446004, -0.8416212335729143, -0.5244005127080409, -0.2533471031357997, 0.0, 0.2533471031357997, 0.5244005127080407, 0.8416212335729143, 1.2815515655446004]
[0.019970984035859413, 0.05805826175840775, 0.1869504831500295, 0.5, 0.8130495168499705, 0.9419417382415922, 0.9800290159641406]
[0.0, 0.0, 0.0, 0.5, 0.3032653298563167, 0.18393972058572114, 0.11156508007421491]
Out[6]:
Normal{Float64}(μ=-0.08351240461430262, σ=1.3479498728939758)
In [7]:
# Statistical Tests in Julia
# Listing J.7
# Last updated July 2020
#
#

using Random, HypothesisTests; # loading required packages

Random.seed!(100)                # set random seed        
x = rand(TDist(5), 500)          # create hypothetical dataset x

println(JarqueBeraTest(x))       # Jarque-Bera test for normality
println(LjungBoxTest(x,20))      # Ljung-Box test for serial correlation
Jarque-Bera normality test
--------------------------
Population details:
    parameter of interest:   skewness and kurtosis
    value under h_0:         "0 and 3"
    point estimate:          "-0.02625371325018951 and 5.7924335679490655"

Test summary:
    outcome with 95% confidence: reject h_0
    one-sided p-value:           <1e-35

Details:
    number of observations:         500
    JB statistic:                   162.509

Ljung-Box autocorrelation test
------------------------------
Population details:
    parameter of interest:   autocorrelations up to lag k
    value under h_0:         "all zero"
    point estimate:          NaN

Test summary:
    outcome with 95% confidence: fail to reject h_0
    one-sided p-value:           0.9101

Details:
    number of observations:         500
    number of lags:                 20
    degrees of freedom correction:  0
    Q statistic:                    12.171

In [8]:
# Time Series in Julia
# Listing J.8
# Last updated July 2020
#
#

using Plots;

Random.seed!(100)
x = rand(TDist(5), 60)           # Create hypothetical dataset x

acf = autocor(x, 1:20)           # autocorrelation for lags 1:20
pacf = autocor(x, 1:20)          # partial autocorrelation for lags 1:20

## Plotting using Plots.jl
plot(bar(acf, title = "Autocorrelation", legend = false), bar(pacf, title = "Partial autocorrelation", legend = false))
Out[8]:
In [9]:
# Loops and Functions in Julia
# Listing J.9
# Last updated July 2020
#
#


## For loops

for i in range(3,length = 5)        # using range with the "length" option
    println(i^2)                    # where n = number of terms      
    end                             # this iterates over [3,4,5,6,7]

## If-else loops

X = 10

if X % 3 == 0
    println("X is a multiple of 3")
else
    println("X is not a multiple of 3")
end

## Functions (example: a simple excess kurtosis function)

using Statistics;

function excess_kurtosis(x, excess = 3)::Float64      # excess optional, default = 3
    m4 = mean((x .- mean(x)).^4)                      # element-wise exponentiation .^
    excess_kurt = m4/(std(x)^4) - excess
    return excess_kurt
end

using Random, Distributions;

Random.seed!(100)
x = rand(TDist(5), 60)           # Create hypothetical dataset x

excess_kurtosis(x)               

## Note: we have forced output to be of type Float64 by the type declaration above
9
16
25
36
49
X is not a multiple of 3
Out[9]:
4.753425237265262
In [10]:
# Basic Graphs in Julia
# Listing J.10
# Last updated July 2020
#
#

## For the simple plots in FRF we use Plots.jl 
## Full documentation at http://docs.juliaplots.org/latest/

## By default, Plots.jl uses the GR backend, sufficient for plots done in FRF
## Alternative backends are also available, e.g. Plotly, PlotlyJS

y = rand(Normal(0,1), 50)

## Plotting barplot, lineplot, histogram, scatterplot of y

return plot(bar(y, title = "Bar plot"), plot(y, title = "Line plot"), 
    histogram(y, title = "Histogram"), scatter(y, title = "Scatter plot"), legend = false)

## Wrapping plot(...) around multiple plots allows for automatic subplotting
## Options in wrapped plot(...) apply to all subplots
## Plot individual graphs using histogram(y), bar(y), etc. directly

## More examples using GR (plus syntax for customizations) can be found online:
## https://docs.juliaplots.org/latest/generated/gr/
Out[10]:
In [11]:
# Miscellaneous Useful Functions in Julia
# Listing J.11
# Last updated July 2020
#
#

## 1) To convert objects from one type to another, use convert(Type, object)
##    To check type, use typeof(object)

x = 8.0

println(typeof(x))

x = convert(Int, 8.0)

println(typeof(x))

## 2) To type Greek letters, type \ + name + Tab in succession
##    e.g. \gammaTab gives you γ and \GammaTab gives you Γ 
##
##    Greek letters are sometimes essential in retrieving parameters from functions
##    e.g. res = fit_mle(Normal, x) will return an object res of type Distribution
##    with fitted parameters res.μ and res.σ

y = rand(Normal(0,1), 100)
res = fit_mle(Normal, y)
println("Fitted mean: ", res.μ)
println("Fitted sd: ", res.σ)
Float64
Int64
Fitted mean: 0.044317955654009006
Fitted sd: 1.1058746524691687


Chapter 1: Financial Markets, Prices and Risk

  • 1.1/1.2: Loading hypothetical stock prices, converting to returns, plotting returns
  • 1.3/1.4: Summary statistics for returns timeseries
  • 1.5/1.6: Autocorrelation function (ACF) plots, Ljung-Box test
  • 1.7/1.8: Quantile-Quantile (QQ) plots
  • 1.9/1.10: Correlation matrix between different stocks
In [12]:
# Download S&P 500 data in Julia
# Listing 1.1/1.2
# Last updated June 2022
#
#

using CSV, DataFrames;

price = CSV.read("index.csv",DataFrame)
y = diff(log.(price[:,1]))

using Plots;
plot(y, title = "S&P500 returns", legend = false)
Out[12]:
In [13]:
# Sample statistics in Julia
# Listing 1.3/1.4
# Last updated July 2020
#
#

using Statistics, StatsBase;

println("Standard deviation: ", std(y), "\n")
println("Minimum value: ", minimum(y), "\n")
println("Maximum value: ", maximum(y), "\n")
println("Skewness: ", skewness(y), "\n")
println("Kurtosis: ", kurtosis(y), "\n")
println("Autocorrelation of returns:", "\n", autocor(y, 1:20), "\n")
println("Autocorrelation of returns squared:", "\n", autocor(y.^2, 1:20), "\n")

using HypothesisTests;

println(JarqueBeraTest(y))
println(LjungBoxTest(y,20))
println(LjungBoxTest(y.^2, 20))
Standard deviation: 0.010005728643763064

Minimum value: -0.10195548627302298

Maximum value: 0.10673589502911707

Skewness: 0.15263326989633666

Kurtosis: 13.981171461092032

Autocorrelation of returns:
[0.01481545997603839, -0.006631341774916777, 0.00897934207314229, 0.038936627359256126, -0.04202186609969318, 0.023510498020909366, 0.04752702011094282, 0.01130212301720032, -0.03975424081522125, -0.012375328429854171, 0.01822253865468312, 0.015346642354154541, -0.026040713086649975, -0.007244542847021728, 0.03619442782751118, -0.04581690373803026, -0.0071677689174158175, 0.04898928881836514, 0.014602826246358232, -0.027598158783722564]

Autocorrelation of returns squared:
[0.19918747286311722, 0.14460509432224403, 0.12364476033877878, 0.20469155789639604, 0.16712735240957596, 0.193404279042487, 0.0791799524328806, 0.11944059864788954, 0.05901632918874656, 0.15627253944262184, 0.13966081960880813, 0.14301595584997348, 0.1191919007014663, 0.12036386334601584, 0.11593141643561286, 0.13805754069254889, 0.2285290987347057, 0.14272081972084968, 0.11620336638371424, 0.16816172360597395]

Jarque-Bera normality test
--------------------------
Population details:
    parameter of interest:   skewness and kurtosis
    value under h_0:         "0 and 3"
    point estimate:          "0.15263326989633633 and 16.98117146109203"

Test summary:
    outcome with 95% confidence: reject h_0
    one-sided p-value:           <1e-99

Details:
    number of observations:         5676
    JB statistic:                   46251.4

Ljung-Box autocorrelation test
------------------------------
Population details:
    parameter of interest:   autocorrelations up to lag k
    value under h_0:         "all zero"
    point estimate:          NaN

Test summary:
    outcome with 95% confidence: reject h_0
    one-sided p-value:           <1e-10

Details:
    number of observations:         5676
    number of lags:                 20
    degrees of freedom correction:  0
    Q statistic:                    93.4881

Ljung-Box autocorrelation test
------------------------------
Population details:
    parameter of interest:   autocorrelations up to lag k
    value under h_0:         "all zero"
    point estimate:          NaN

Test summary:
    outcome with 95% confidence: reject h_0
    one-sided p-value:           <1e-99

Details:
    number of observations:         5676
    number of lags:                 20
    degrees of freedom correction:  0
    Q statistic:                    2543.05

In [14]:
# ACF plots and the Ljung-Box test in Julia
# Listing 1.5/1.6
# Last updated July 2020
#
#


q1 = autocor(y, 1:20)
q2 = autocor(y.^2, 1:20)

plot(bar(q1, title = "ACF of returns"), 
    bar(q2, title = "ACF of returns squared"), legend = false)
Out[14]:
In [15]:
# QQ plots in Julia
# Listing 1.7/1.8
# Last updated July 2020
#
#

using StatsPlots, Distributions;

plot(qqplot(Normal,float(y),qqline=:quantile, title = "QQPlot vs Normal"), 
    qqplot(TDist(5),float(y),qqline=:quantile, title = "QQPlot vs Student-t(5)"),
    qqplot(TDist(4),float(y),qqline=:quantile, title = "QQPlot vs Student-t(4)"),    
    qqplot(TDist(3),float(y),qqline=:quantile, title = "QQPlot vs Student-t(3)"))
Out[15]: