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

#--------------------------------------------------------------------------
# Include procedures for model solution
#--------------------------------------------------------------------------
readDirSol = "$(pwd())/FunctionsSolution/"
include(readDirSol * "cga_grid_h.jl")
include(readDirSol * "cga_grid_km.jl")
include(readDirSol * "cga_grid_kms.jl")
include(readDirSol * "eval_hh_33_revised.jl")
include(readDirSol * "get_coefs_plc.jl")
include(readDirSol * "get_GH_nodes.jl")
include(readDirSol * "get_order1_sgu.jl")
include(readDirSol * "ghpoints.jl")
include(readDirSol * "gx_hx.jl")
include(readDirSol * "my_levenberg_marquardt2.jl")
include(readDirSol * "options.jl")
include(readDirSol * "PLC_solve.jl")
include(readDirSol * "PLC_main.jl")
#include(readDirSol * "simulate_plc.jl")
include(readDirSol * "solve_sgu.jl")
#include(readDirSol * "system_piecewise_linear_jacob_33_revised.jl")
include(readDirSol * "tsga_grid.jl")
include(readDirSol * "Monomials_2_N4.jl")

#--------------------------------------------------------------------------
# Include procedures for model estimation
#--------------------------------------------------------------------------
readDir = "FunctionsEstimation/"
include(readDir * "fCreateParamDict.jl")
include(readDir * "fAddSteadyStateParamDict.jl")
include(readDir * "fCanonicalForm.jl")
include(readDir * "fStateSpaceForm.jl")
include(readDir * "fStateSpaceSimulationaug.jl")
include(readDir * "fCOPFaugExactZLB.jl")
include(readDir * "fCondDistExactZLB_v5.jl")
include(readDir * "fPriorDraws.jl")
include(readDir * "fPriorDensity.jl")
include(readDir * "run_filter_COPFexactZLB.jl")
include(readDir * "fSystematicResampling.jl")

nM        = 180; # number of particles
ME_scale  = 0.001;;
mhrun     = 1;
silent    = 1;
nSistResFilter = 0; # = 1 systematic resampling in filter, multinomial otherwise
pZLB      = 0;
rnd_off   = 0;

#--------------------------------------------------------------------------
# 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
   mData = convert(Array, mData);
   R_quantile = quantile(mData[nStartDate:nEndDate,5]./400 .+ 1,0.90) # R data quantile

elseif  nDataSet == "US2007Q4"
   mData    = CSV.read("Data/ACS-RED-USDATA-G.csv")
   nStartDate = 1
   nEndDate   = 96 # 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

#--------------------------------------------------------------------------
# Set solution type and Grid
#--------------------------------------------------------------------------
nSolutionRCLinear = 1
nSmolyakGrid = 1 # = 1 use Smolyak Grid. Note: need to set in options.jl file options_int_out["ext"] = 1 and options_int_out["k"] = 0

options_int_use,options_float_use = @eval options()

if nSmolyakGrid == 1
    options_int_use["k"]   = 0 # Number of non-linear solve steps
    options_int_use["ext"] = 1 # 1 = load cluster from input

    include("FunctionsSolutionSmolyak/Smolyak_Elem_isotrop.jl")
    include("FunctionsSolutionSmolyak/Smolyak_Grid.jl")
    include("FunctionsSolutionSmolyak/fGenSmolyakGrid.jl")
    include("FunctionsSolutionSmolyak/fScaleSmolyakGrid2.jl")

    d          = 6
    mu         = 2
    sd         = 1.282
    wSmolyak   = 1

    vSmolyakGridSpec = [d, mu, sd, wSmolyak, R_quantile]
    fixed_grid = fGenSmolyakGrid(d,mu,sd)

else
    fixed_grid = zeros(1,1)
    vSmolyakGridSpec = ones(1)
end

if nSolutionRCLinear == 1
    include(readDirSol * "simulate_plc_lin.jl")
    include(readDirSol * "system_piecewise_linear_jacob_33_lin.jl")
else
    include(readDirSol * "simulate_plc_nlin.jl")
    include(readDirSol * "system_piecewise_linear_jacob_33_nlin.jl")
end

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

nFileStrMH   = "COPFexactZLB_" * nNamePriorSpec * "_" * nDataSet * "_RClin" * "$(nSolutionRCLinear)" * nGrid * "_Mhrun" * "$(mhrun)"
nFileStrFilt = "COPFexactZLB_Filtered_" * nNamePriorSpec * "_" * nDataSet * "_RClin" * "$(nSolutionRCLinear)" * nGrid * "_Mhrun" * "$(mhrun)"
#nFileStrFilt =  nDataSet * "_RClin" * "$(nSolutionRCLinear)" * nGrid * "_Mhrun" * "$(mhrun)"
savedir = "$(pwd())/PosteriorDraws/" * nFileStrMH * "/Filtered/"
try mkdir(savedir) catch; end  # Create directory, if it exists skip creation

run_filter_COPFexactZLB(ME_scale, nM, mhrun, silent, mData, nFileStrMH, nFileStrFilt, nSolutionRCLinear,
                        rnd_off,fixed_grid,vSmolyakGridSpec,options_int_use,options_float_use)
