 # Appendix - Introduction (in MATLAB/Julia)

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/.

##### Listing M.1: Entering and Printing Data Last updated June 2018

x = 10; % assign x the value 10, silencing output print with ;
disp(x) % display x

##### Listing J.1: Entering and Printing Data in Julia 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


##### Listing M.2: Vectors, Matrices and Sequences Last updated June 2018

y = [1,3,5,7,9]      % lists are denoted by square brackets
y(3)       % calling 3rd element (MATLAB indices start at 1)
size(y)    % shows that y is 1 x 5 (a row vector, by default)
length(y)  % as expected, y has length 5
v = nan(2,3)         % fill a 2 x 3 matrix with NaN values
size(v)    % as expected, v is size (2,3)
w = repmat([1,2,3]', 2, 3) % repeats matrix twice by rows, thrice by columns
s = 1:10   % s is a list of integers from 1 to 10 inclusive

##### Listing J.2: Vectors, Matrices and Sequences in Julia Last updated July 2020

y = [1,3,5,7,9]      # lists in square brackets are stored as arrays
println(y)
println(y)        # 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


##### Listing M.3: Basic Summary Statistics Last updated June 2022

y = [3.14; 15; 9.26; 5]; % List with semicolons is a column vector
sum(y)     % sum of all elements of y
prod(y)    % product of all elements of y
max(y)     % maximum value of y
min(y)     % minimum value of y
range(y)   % difference between max and min value of y
mean(y)    % arithmetic mean
median(y)  % median
var(y)     % variance
cov(y)     % covar matrix = variance for single vector
corrcoef(y)          % corr matrix =  for single vector
sort(y)    % sorting in ascending order
log(y)     % natural log

##### Listing J.3: Basic Summary Statistics in Julia 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 =  for single vector
println(sort(y))     # sorts y in ascending order
println(log.(y))     # natural log, note . denotes elementwise operation


##### Listing M.4: Calculating Moments Last updated June 2018

mean(y)     % mean
var(y)     % variance
std(y)     % unbiased standard deviation, by default
skewness(y) % skewness
kurtosis(y) % kurtosis

##### Listing J.4: Calculating Moments in Julia 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)


##### Listing M.5: Basic Matrix Operations Last updated June 2018

z = [1, 2; 3, 4] % z is a 2 x 2 matrix (Note the use of ; as row separator)
x = [1, 2]           % x is a 1 x 2 matrix
%% Note: z * x is undefined since the two matrices are not conformable
z * x'     % this evaluates to a 2 x 1 matrix
vertcat(z,x)     % "stacking" z and x vertically
horzcat(z,x')    % "stacking z and x' horizontally
%% Note: dimensions must match along the combining axis)

##### Listing J.5: Basic Matrix Operations in Julia 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


##### Listing M.6: Statistical Distributions Last updated June 2018

q = -3:1:3           % specify a set of values
p = 0.1:0.1:0.9      % specify a set of probabilities
norminv(p, 0, 1)     % element-wise inverse Normal quantile
tcdf(q, 4)           % element-wise cdf under Student-t(4)
chi2pdf(q, 2)        % element-wise pdf under Chisq(2)
%% One can also obtain pseudorandom samples from distributions
x = trnd(5, 100, 1);           % Sampling 100 times from t dist with 5 df
y = normrnd(0, 1, 100, 1); % Sampling 50 times from a standard normal
%% Given sample data, we can also obtain MLE estimates of distribution parameters:
res = fitdist(x, "Normal") % Fitting x to normal dist

##### Listing J.6: Statistical Distributions in Julia 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/


##### Listing M.7: Statistical Tests Last updated July 2020

x = trnd(5, 500, 1);           % Create hypothetical dataset x
[h1, p1, jbstat] = jbtest(x)   % Jarque-Bera test for normality
[h2, p2, lbstat] = lbqtest(x,'lags',20) % Ljung-Box test for serial correlation - Needs Econometrics Toolbox

##### Listing J.7: Statistical Tests in Julia Last updated July 2020

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


##### Listing M.8: Time Series Last updated June 2018

x = trnd(5, 60, 1); % Create hypothetical dataset x
subplot(1,2,1)
autocorr(x, 20)     % autocorrelation for lags 1:20
subplot(1,2,2)
parcorr(x,20)        % partial autocorrelation for lags 1:20

##### Listing J.8: Time Series in Julia 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))


##### Listing M.9: Loops and Functions Last updated June 2018

%% For loops
for i = 3:7 % iterates through [3,4,5,6,7]
i^2
end
%% If-else loops
X = 10;
if (rem(X,3) == 0)
disp("X is a multiple of 3")
else
disp("X is not a multiple of 3")
end
%% Functions (example: a simple excess kurtosis function)
%% NOTE: in MATLAB, functions can be defined in 2 locations:
%% 1) in a separate file (e.g. excess_kurtosis.m in this case) in the workspace
%% 2) in the same file as the rest of the code, BUT at the end of the file
%% function k = excess_kurtosis(x, excess)
%%     if nargin == 1 % if there is only 1 argument
%%         excess = 3; % set excess = 3
%%     end  % this is how optional param excess is set
%%     m4 = mean((x-mean(x)).^4);
%%     k = m4/(std(x)^4) - excess;
%% end

##### Listing J.9: Loops and Functions in Julia 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


##### Listing M.10: Basic Graphs Last updated July 2020

y = normrnd(0, 1, 50, 1);
z = trnd(4, 50, 1);
subplot(2,2,1)
bar(y)     % bar plot
title("Bar plot")
subplot(2,2,2)
plot(y)    % line plot
title("Line plot")
subplot(2,2,3)
histogram(y) % histogram
title("Histogram")
subplot(2,2,4)
scatter(y,z) % scatter plot
title("Scatter plot")

##### Listing J.10: Basic Graphs in Julia 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/


##### Listing M.11: Miscellaneous Useful Functions Last updated June 2018

%% Convert objects from one type to another with int8() etc
%% To check type, use isfloat(object), isinteger(object) and so on
x = 8.0;
isfloat(x)
x = int8(x);
isinteger(x)

##### Listing J.11: Miscellaneous Useful Functions in Julia 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.σ)