R and Julia Appendix - Introduction

Appendix - Introduction

R and Julia

Copyright 2011 - 2023 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: www.gnu.org/licenses.

Listing 0.1
Entering and Printing Data in R
x = 10             # assign x the value 10
print(x)           # print x
Listing 0.1
Entering and Printing Data in Julia
x = 10             # assign x the value 10
println(x)         # print x

Listing 0.2
Vectors, Matrices and Sequences in R
y = c(1,3,5,7,9)          # create vector using c()
print(y)
print(y[3])               # calling 3rd element (R indices start at 1)
print(dim(y))             # gives NULL since y is a vector, not a matrix
print(length(y))          # as expected, y has length 5
v = matrix(nrow=2,ncol=3) # fill a 2 x 3 matrix with NaN values (default)
print(dim(v))             # as expected, v is size (2,3)
w = matrix(c(1,2,3),nrow=6,ncol=3) # repeats matrix twice by rows, thrice by columns
print(w)
s = 1:10                  # s is a list of integers from 1 to 10 inclusive
print(s)                  
Listing 0.2
Vectors, Matrices and Sequences in Julia
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

Listing 0.3
Basic Summary Statistics in R
y=matrix(c(3.1,4.15,9))
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)       # min, max value of y
mean(y)        # arithmetic mean
median(y)      # median
var(y)         # variance
cov(y)         # covar matrix = variance for single vector
cor(y)         # corr matrix = [1] for single vector
sort(y)        # sorting in ascending order
log(y)         # natural log
Listing 0.3
Basic Summary Statistics in Julia
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

Listing 0.4
Calculating Moments in R
library(moments)
mean(y)      # mean
var(y)       # variance
sd(y)        # unbiased standard deviation, by default
skewness(y)  # skewness
kurtosis(y)  # kurtosis
Listing 0.4
Calculating Moments in Julia
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 0.5
Basic Matrix Operations in R
z = matrix(c(1,2,3,4),2,2)   # z is a 2 x 2 matrix
x = matrix(c(1,2),1,2)       # x is a 1 x 2 matrix
z %*% t(x)                   # this evaluates to a 2 x 1 matrix
rbind(z,x)                   # "stacking" z and x vertically
cbind(z,t(x))                # "stacking z and x' horizontally
Listing 0.5
Basic Matrix Operations in Julia
z = Matrix([[1 2];[3 4]])   # z is a 2 x 2 matrix
x = Matrix([1 2])           # x is a 1 x 2 matrix
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 

Listing 0.6
Statistical Distributions in R
q = seq(from = -3, to = 3, length = 7)     # specify a set of values
p = seq(from = 0.1, to = 0.9, length = 9)  # specify a set of probabilities
qnorm(p, mean = 0, sd = 1)                 # element-wise inverse Normal quantile
pt(q, df = 4)                              # element-wise cdf under Student-t(4)
dchisq(q, df = 2)                          # element-wise pdf under Chisq(2)
x = rt(100, df = 5)                        # Sampling 100 times from TDist with 5 df
y = rnorm(50, mean = 0, sd = 1)            # Sampling 50 times from a standard normal 
library(MASS)
res = fitdistr(x, densfun = "normal")      # Fitting x to normal dist
print(res)
Listing 0.6
Statistical Distributions in Julia
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)
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 
fit_mle(Normal, x)                 # Fitting x to normal dist     

Listing 0.7
Statistical Tests in R
library(tseries)
x = rt(500, df = 5)                            # Create hypothetical dataset x
jarque.bera.test(x)                            # Jarque-Bera test for normality
Box.test(x, lag = 20, type = c("Ljung-Box"))   # Ljung-Box test for serial correlation
Listing 0.7
Statistical Tests in Julia
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

Listing 0.8
Time Series in R
x = rt(60, df = 5)  # Create hypothetical dataset x
par(mfrow=c(1,2), pty='s')
acf(x,20)           # autocorrelation for lags 1:20
pacf(x,20)          # partial autocorrelation for lags 1:20
Listing 0.8
Time Series in Julia
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
plot(bar(acf, title = "Autocorrelation", legend = false))
plot(bar(pacf, title = "Partial autocorrelation", legend = false)) 

Listing 0.9
Loops and Functions in R
for (i in 3:7)        # iterates through [3,4,5,6,7]
    print(i^2)      
X = 10
if (X %% 3 == 0) {
    print("X is a multiple of 3")
} else {
    print("X is not a multiple of 3")
}
excess_kurtosis = function(x, excess = 3){ # note: excess optional, default=3
    m4 = mean((x-mean(x))^4)
    excess_kurt = m4/(sd(x)^4) - excess
    excess_kurt
}
x = rt(60, df = 5)                         # Create hypothetical dataset x
excess_kurtosis(x)      
Listing 0.9
Loops and Functions in Julia
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]
X = 10
if X % 3 == 0
    println("X is a multiple of 3")
else
    println("X is not a multiple of 3")
end
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)               

Listing 0.10
Basic Graphs in R
y = rnorm(50, mean = 0, sd = 1)
par(mfrow=c(2,2)) # sets up space for subplots
barplot(y)        # bar plot
plot(y,type='l')  # line plot
hist(y)           # histogram
plot(y)           # scatter plot
Listing 0.10
Basic Graphs in Julia
y = rand(Normal(0,1), 50)
return plot(bar(y, title = "Bar plot"), plot(y, title = "Line plot"), 
    histogram(y, title = "Histogram"), scatter(y, title = "Scatter plot"), legend = false)

Listing 0.11
Miscellaneous Useful Functions in R
x = 8.0
print(typeof(x))
x = as.integer(x)
print(typeof(x))
Listing 0.11
Miscellaneous Useful Functions in Julia
x = 8.0
println(typeof(x))
x = convert(Int, 8.0)
println(typeof(x))
y = rand(Normal(0,1), 100)
res = fit_mle(Normal, y)
println("Fitted mean: ", res.μ)
println("Fitted sd: ", res.σ)


Financial Risk Forecasting
Market risk forecasting with R, Julia, Python and Matlab. Code, lecture slides, implementation notes, seminar assignments and questions.
© All rights reserved, Jon Danielsson,