Appendix - Introduction

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

y = [1,3,5,7,9]       # lists in square brackets are stored as arrays
print(y)
print(y[2])           # 3rd element (Python indices start at 0)
print(len(y))         # as expected, y has length 5
import numpy as np    # NumPy: Numeric Python package
v = np.full([2,3], np.nan) # create a 2x3 matrix with NaN values
print(v)           
print(v.shape)        # as expected, v is size (2,3)
w = np.tile([1,2,3], (3,2)) # repeats thrice by rows, twice by columns
print(w)
s = range(10)         # an iterator from 0 to 9
print([x for x in s]) # return  elements using list comprehension

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

import numpy as np
y = [3.14, 15, 9.26, 5]
print(sum(y))          # sum of all elements of y
print(max(y))          # maximum value of y
print(min(y))          # minimum value of y
print(np.mean(y))      # arithmetic mean
print(np.median(y))    # median
print(np.var(y))       # population variance
print(np.cov(y))       # covar matrix = sample variance for single vector
print(np.corrcoef(y))  # corr matrix = [1] for single vector
print(np.sort(y))      # sort in ascending order
print(np.log(y))       # natural log

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

from scipy import stats
print(np.mean(y))                        # mean
print(np.var(y))                         # variance
print(np.std(y, ddof = 1))               # ddof = 1 for unbiased standard deviation
print(stats.skew(y))                     # skewness
print(stats.kurtosis(y, fisher = False)) # fisher = False gives Pearson definition

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)

z = np.matrix([[1, 2], [3, 4]])                   # z is a 2 x 2 matrix
x = np.matrix([1, 2])                             # x is a 1 x 2 matrix
print(z * np.transpose(x))                        # this evaluates to a 2 x 1 matrix          
b = np.concatenate((z,x), axis = 0)               # "stacking" z and x vertically
print(b)
c = np.concatenate((z,np.transpose(x)), axis = 1) # "stacking" z and x horizontally
print(c)

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 

q = np.arange(-3,4,1)        # specify a set of values, syntax arange(min, exclusive-max, step)
p = np.arange(0.1,1.0,0.1)   # specify a set of probabilities
print(stats.norm.ppf(p))     # element-wise inverse Normal quantile
print(stats.t.cdf(q,4))      # element-wise cdf under Student-t(4)
print(stats.chi2.pdf(q,2))   # element-wise pdf under Chisq(2)
x = np.random.standard_t(df=5, size=100)   # Sampling 100 times from TDist with 5 df
y = np.random.normal(size=50)              # Sampling 50 times from a standard normal 
res = stats.norm.fit(x)                    # Fitting x to normal dist
print(res)                                 # First element is mean, second sd

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     

from statsmodels.stats.diagnostic import acorr_ljungbox
x = np.random.standard_t(df=5, size=500)  # Create dataset x
print(stats.jarque_bera(x))               # Jarque-Bera test - prints statistic and p-value
print(acorr_ljungbox(x, lags=20))         # Ljung-Box test - prints array of statistics and p-values

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

import statsmodels.api as sm
import matplotlib.pyplot as plt
y = np.random.standard_t(df = 5, size = 60)   # Create hypothetical dataset y
q1 = sm.tsa.stattools.acf(y, nlags=20)        # autocorrelation for lags 1:20
plt.bar(x = np.arange(1,len(q1)), height = q1[1:])
plt.show()
plt.close()
q2 = sm.tsa.stattools.pacf(y, nlags=20)       # partial autocorr for lags 1:20
plt.bar(x = np.arange(1,len(q2)), height = q2[1:])
plt.show()
plt.close()

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)) 

for i in range(3,8):      # NOTE: range(start, end), end excluded
    print(i**2)           # range(3,8) iterates through [3,4,5,6,7)
X = 10
if X % 3 == 0:
    print("X is a multiple of 3")
else:
    print("X is not a multiple of 3")
def excess_kurtosis(x, excess = 3):        # note: excess optional, default = 3
    m4=np.mean((x-np.mean(x))**4)          # note: exponentiation in Python uses **
    excess_kurt=m4/(np.std(x)**4)-excess
    return excess_kurt
x = np.random.standard_t(df=5,size=60)     # Create hypothetical dataset x
print(excess_kurtosis(x))

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)               

y = np.random.normal(size = 50)
z = np.random.standard_t(df = 4, size = 50)
plt.subplot(2,2,1)
plt.bar(range(len(y)), y);
plt.subplot(2,2,2)
plt.plot(y);
plt.subplot(2,2,3)
plt.hist(y);
plt.subplot(2,2,4)
plt.scatter(y,z);

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)

x = 8.0
print(type(x))
x = int(x)
print(type(x))

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.σ)