
using Atom
using Distributions
using Clustering
using Optim
#using LsqFit
using Distances
using StatsBase
using DataFrames
using CSV
using Random
using LinearAlgebra
#using ControlSystems # for discrete Lyapunov equation solver

clearconsole()

#--------------------------------------------------------------------------
# Include procedures for model solution
#--------------------------------------------------------------------------
readDir = "$(pwd())/FunctionsSolution/"
include(readDir * "Linear_main.jl")
include(readDir * "solve_sgu.jl")
include(readDir * "options.jl")
include(readDir * "get_order1_sgu.jl")
include(readDir * "gx_hx.jl")

#--------------------------------------------------------------------------
# Include procedures for model estimation
#--------------------------------------------------------------------------
readDir = "$(pwd())/FunctionsEstimation/"
include(readDir * "fCreateParamDict.jl")
include(readDir * "fAddSteadyStateParamDict.jl")
include(readDir * "fCanonicalFormLinear.jl")
include(readDir * "fStateSpaceFormLinear.jl")
include(readDir * "fStateSpaceSimulationLinear.jl")
include(readDir * "fPriorDensity.jl")
include(readDir * "fPriorDraws.jl")
include(readDir * "run_MCMC_KF.jl")
include(readDir * "fKF.jl")
include(readDir * "fCanonicalFormLinear.jl")

ME_scale  = 0.0;
nSimRWMH  = 110000;
mhrun     = 1;
nBurnRWMH = 10000;      # number of simulations for the Metropolis-Hasting algorithm
silent    = 1;
nOutput   = 10000;      # print output every ncount steps
c         = 0.3;     # scaling of proposal covariance matrix for Sim1 data
                       # use 0.002 if constructed from Prior
                       # use 0.3 if constructed from Posterior
nSistResFilter = 0; # = 1 systematic resampling in filter, multinomial otherwise

#--------------------------------------------------------------------------
# Choose prior
#--------------------------------------------------------------------------

nNamePriorSpec = "Prior1"

#--------------------------------------------------------------------------
# Load data set
# date obs ygrowth (%) infl (annual %) R (annual %)
# 100*logcy_raw 100*logcy 100*logcy4
#--------------------------------------------------------------------------
nDataSet = "US" # US or Sim1

if nDataSet == "US"
   mData    = CSV.read("Data/ACS-RED-USDATA-G.csv")
   nStartDate = 1
   nEndDate   = 140 # 140 is full sample, 96 is preZLB
   mData = convert(Array, mData);
   R_quantile = quantile(mData[nStartDate:nEndDate,5]./400 .+ 1,0.90) # R data quantile

else # if nDataSet == "Sim1"
   mData    = CSV.read("Data/COPFexactZLB_DGP1_Data.csv")
   nStartDate = 401
   nEndDate   = 540
   mData = convert(Array, mData);
   R_quantile = 1.01

end

nFileStr   = "KF_" * nNamePriorSpec * "_" * nDataSet * "_RClin1_GridNone"
nFileStrMH = "KF_" * nNamePriorSpec * "_" * nDataSet * "_RClin1_GridNone" * "_Mhrun" *"$(mhrun)"
savedir      = "$(pwd())/PosteriorDraws/" * nFileStrMH * "/"
try mkdir(savedir) catch; end  # Create directory, if it exists skip creation
#--------------------------------------------------------------------------
# Construct proposal covariance matrix
#--------------------------------------------------------------------------

nLoadCovVar = 1;

if nLoadCovVar == 0
    readDir       = "$(pwd())/PriorDraws/"
    readName      = "SummaryDirectDrawsPrior1.csv"
    mSummaryDraws = CSV.read(readDir * readName ; datarow=1)
    mSummaryDraws = convert(Array,mSummaryDraws)
    vStdDevTheta  = mSummaryDraws[2,:]
    mSigmaProp    = Diagonal(vStdDevTheta.^2)
else
    readDir       = "$(pwd())/PosteriorDraws/" * nFileStr * "_Mhrun0/"
    readName      = nFileStr * "_Mhrun0_CovVarMatrix.csv"
    mSigmaProp    = CSV.read(readDir * readName ; datarow=1)
    mSigmaProp    = convert(Array,mSigmaProp)
end

#--------------------------------------------------------------------------
# Initial Value for MCMC
#--------------------------------------------------------------------------

readDir      = "$(pwd())/PosteriorDraws/"
readName     = nFileStrMH * "_InitialValue.csv"
vParameters  = CSV.read(readDir * readName , datarow=1 , limit=1)
vParameters  = convert(Array,vParameters)'


run_MCMC_KF(ME_scale, nSimRWMH, mhrun, nBurnRWMH,
           silent, nOutput, c, nNamePriorSpec, mData, nStartDate, nEndDate,
           nFileStr, nFileStrMH, vParameters, mSigmaProp)
