Multivariate-Normal-Dist-Examples-Julia.jl

## MIT license (C) 2019 by Andrew Lyasoff

# Julia 1.1.0 code to illustrate the construction of univariate histograms
# and the Monte Carlo simulation technique for multivariate Gaussian samples with a given covariance matrix.
#See (3.67), p.91, in "Stochastic Methods in Asset Pricing."

begin
    using LinearAlgebra
    using SpecialFunctions
    using StatsBase
    using Random
    using Plots
    pyplot()
end

#We write our own histogram function. It takes as an input a single 1-dimensional array of data.
#The number of bins in the histogram is determined automatically by using the Diaconis-Friedman rule.
#The function returns two arrays: the mid-points of the bins and the (unnormalized) heights of the bars.

function hstgram(data_sample::Array{Float64,1})
    data_sorted=sort(data_sample)
    first=data_sorted[1]
    last=data_sorted[end]
    nmb=length(data_sorted)
    IQR=percentile(data_sorted,75)-percentile(data_sorted,25)
    bin_size_loc = 2*IQR*(nmb^(-1.0/3))
    num_bins=Int(floor((last-first)/bin_size_loc))
    bin_size=(last-first)/(num_bins)
    bin_end_points=[first+(i-1)*bin_size for i=1:(num_bins+1)]
    ahist_val=[length(data_sorted[data_sorted .< u]) for u in bin_end_points]
    hist_val=[ahist_val[i+1]-ahist_val[i] for i=1:num_bins]
    mid_bins=[first-bin_size/2+i*bin_size for i=1:num_bins]
    return mid_bins, hist_val
end


## First, create data sampled from the standard univariate normal density.
val=(x->((2*π)^(-1/2)*exp(-x^2/2))).(-3.3:0.05:3.3);

Random.seed!(0xabcdef12); # if needed to generate the same samples

#simulate 10,000 standard normals and generate the histogram
begin
    nval=randn!(zeros(10000));
    U,V=hstgram(nval);
    VV=V/(sum(V)*(U[2]-U[1]));
end

#check the length of the bin
length(VV)

begin
    plot(U.+(U[2]-U[1])/2,VV,line=(:steppre,1),linewidth=0.05,label="histogram")
    xlabel!("samples")
    ylabel!("frequency")
    plot!(-3.3:0.05:3.3,val,label="normal density")
end

#generate another sample
begin
    nval=randn(10000);
    U,V=hstgram(nval);
    VV=V/(sum(V)*(U[2]-U[1]));
    plot(U.+(U[2]-U[1])/2,VV,line=(:steppre,1),linewidth=0.05,label="histogram")
    xlabel!("samples")
    ylabel!("frequency")
    plot!(-3.3:0.05:3.3,val,label="normal density")
end


# standard normals can be generated by sampling from the uniform distribution
begin
    uval=rand(10000);
    nval=(x->sqrt(2)*erfinv(2*x-1)).(uval);
    U,V=hstgram(nval);
    VV=V/(sum(V)*(U[2]-U[1]));
    plot(U.+(U[2]-U[1])/2,VV,line=(:steppre,1),linewidth=0.05,label="histogram")
    xlabel!("samples")
    ylabel!("frequency")
    plot!(-3.3:0.05:3.3,val,label="normal density")
end

#simulate 10,000 standard bi-variate normals
begin
    nval=randn!(zeros(10000));
    nnval=randn!(zeros(10000));
end

scatter(nval,nnval,ratio=1,markersize=1,label="")

# to simulate bi-variate normals with a given covariance matrix
#   first create a fictitios covariance matrix


begin
    A=rand(2,2);
    #Cov=A'A; # another alternative that always yields a positive definite matrix
    Cov=Symmetric(A) # may not produce a positive definite matrix
end

#check if Cov is positive definite; if not, repeat the last step
eigvals(Cov)

begin
    eigdCov=Diagonal(eigvals(Cov).^0.5)
    eigCov=eigvecs(Cov); # matrix of eigen vectors
    chlCov=cholesky(Cov); # Cholesky "square root" of Cov
    MM=eigdCov*eigCov; # Spectral "square root" of Cov
end

# check that the factorizations give what is expected
Cov-MM'*MM
Cov-chlCov.L*chlCov.U
Cov-(chlCov.U)'*chlCov.U
Cov-(chlCov.L)*(chlCov.L)'

#this should be an orthogonal matrix
eigCov'*eigCov


#Transform the randomly generated standard bi-variate sample through the "square root" of the covariance matrix.
#method 1
begin
    NN=hcat(nval,nnval)';
    data_2_dim=MM'NN;
    scatter(data_2_dim[1,:],data_2_dim[2,:],ratio=1,markersize=1,label="")
end



#method 2
begin
    data_2_dim=(chlCov.L)*NN;
    scatter(data_2_dim[1,:],data_2_dim[2,:],ratio=1,markersize=1,label="")
end