#--------------------------------------------------------------------------
# Script to:
# Generate draws from prior using RWMH
#--------------------------------------------------------------------------

using Atom
using Distributions
using Clustering
using Optim
using LsqFit
using Distances
import StatsBase.sample
import StatsBase.pweights
using CSV
using DataFrames

clearconsole()
#--------------------------------------------------------------------------
# Load the draws
#   each row is a draw for the parameter vector
#   each column is a vector of draws for a specific parameter
#--------------------------------------------------------------------------

nFilter = "BSPF"; # Choose COPFexactZLB or COPF or KF or BSPF
nPrior  = "1";
nMhrun  = "3";
nDataSet = "Sim1" # US Or Sim1
nSolutionRCLinear = 1
nSmolyakGrid = 1


if nSmolyakGrid == 1
    nGrid = "_GridSmolyak"
else
    nGrid = "_GridSim"
end

if nFilter == "KF"
    nGrid = "_GridNone"
    nSolutionRCLinear = 1
end

nFileStrMH = nFilter * "_Prior" * nPrior * "_" * nDataSet * "_RClin" * "$(nSolutionRCLinear)" * nGrid * "_Mhrun" * nMhrun

readDir           = "$(pwd())/PosteriorDraws/" * nFileStrMH * "/"
readNameDraws     = nFileStrMH * "_Draws.csv"
readNameLopPost   = nFileStrMH * "_LogPosterior.csv"

mDraws            = CSV.read(readDir * readNameDraws;header=0)
mDraws = convert(Array,mDraws)

vLogPost          = CSV.read(readDir * readNameLopPost;header=0)
vLogPost = convert(Array,vLogPost)

# Get Max Log Posterior index and maximum posterior density estimator
nMaxLogPOstId = findmax(vLogPost)[2][1]

vMPDEstimator = mDraws[nMaxLogPOstId,:]

# Get CovVar Matrix from draws
mCovMat = cov(mDraws)
saveDir      = readDir
saveName     = nFileStrMH * "_CovVarMatrix.csv"
CSV.write(saveDir * saveName,  DataFrame(mCovMat), writeheader=false)

#--------------------------------------------------------------------------
# Compute means and variances
#--------------------------------------------------------------------------
vMeanTheta   = mean(mDraws; dims = 1)
vStdDevTheta = std(mDraws; dims = 1)

#--------------------------------------------------------------------------
# HPD intervals: take as input posterior draws
#--------------------------------------------------------------------------
include("$(pwd())/FunctionsEstimation/" * "fHPDinterval.jl")
nHPDprob = 0.90;
mHPDinterval = fHPDinterval(mDraws,nHPDprob)

#--------------------------------------------------------------------------
# Save the results
#--------------------------------------------------------------------------
mOutput = [vMeanTheta; vStdDevTheta; mHPDinterval; vMPDEstimator']
saveDir      = readDir
saveName     = nFileStrMH * "_Summary.csv"
CSV.write(saveDir * saveName,  DataFrame(mOutput), writeheader=false)
