
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

clearconsole()

#--------------------------------------------------------------------------
# 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 * "Linear_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 * "fCanonicalFormLinear.jl")
include(readDir * "fStateSpaceForm.jl")
include(readDir * "fStateSpaceFormLinear.jl")
include(readDir * "fStateSpaceSimulationaug.jl")
include(readDir * "fStateSpaceSimulation.jl")
include(readDir * "fBSPF.jl")
include(readDir * "fCOPFaugExactZLB.jl")
include(readDir * "fCOPF.jl")
include(readDir * "fKF.jl")
include(readDir * "fCondDistExactZLB_v5.jl")
include(readDir * "fCondDist_v5.jl")
include(readDir * "fPriorDraws.jl")
include(readDir * "fPriorDensity.jl")
include(readDir * "run_LikelihoodEval.jl")
include(readDir * "fSystematicResampling.jl")

nSimEval   = 100;
nParamDraw = 100;
ME_scale   = 0.05;

nFilter1   = "BSPF" # KF, COPFexactZLB or BSPF
nM_1       = 1000; # number of particles

nFilter2   = "COPF" # KF, COPFexactZLB or BSPF
nM_2       = 180; # number of particles

nFilter3   = "BSPF" # KF, COPFexactZLB or BSPF
nM_3       = 10000; # number of particles

vFilter    = [ nFilter1 nFilter2 ]
vM         = [ nM_1 nM_2 ]


nSistResFilter = 0; # = 1 systematic resampling in filter, multinomial otherwise
pZLB = 1E-6;
rnd_off = 0;

#--------------------------------------------------------------------------
# Load parameter draws to evaluate the filters
#--------------------------------------------------------------------------
#nNamePriorSpec = "Prior1"

mDraw = CSV.read("LikelihoodEvaluations/COPF_Prior1_Sim1_RClin1_GridSmolyak_Mhrun1_SubsampleSpacing200.csv"; header=0)
mDraw = convert(Array,mDraw)
mDraw = mDraw[1:nParamDraw,:]

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

if nDataSet == "Sim1"
   mData    = CSV.read("Data/COPFexactZLB_DGP1_Data.csv")
   mState   = CSV.read("Data/COPFexactZLB_DGP1_States.csv")
   nStartDate = 401
   nEndDate   = 540 # 540
   mData      = convert(Array, mData);
   mState     = convert(Array, mState);
   vInitS     = mState[(nStartDate-1), 3:end]
   R_quantile = 1.01
   vTrueParam = convert(Array,CSV.read("Data/COPFexactZLB_DGP1_Param.csv"; header=0))
   mDraw[1,:] = vTrueParam
   mDraw[:,21:end] = ones(nParamDraw,1) .* vInitS'

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 == "US"
  mData    = CSV.read("Data/ACS-RED-USDATA-G.csv")
  nStartDate = 1
  nEndDate   = 140 # 140 or 95
  mData = convert(Array, mData);
  R_quantile = quantile(mData[nStartDate:nEndDate,5]./400 .+ 1,0.90) # R data quantile

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 # benchmark: 1.282, alt 1.64
    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

run_LikelihoodEval(nSimEval, ME_scale, vFilter, vM, nDataSet,
                mData, nStartDate, nEndDate, mDraw, nParamDraw,
                nSistResFilter,
                pZLB, rnd_off, nSolutionRCLinear, nGrid,
                fixed_grid,vSmolyakGridSpec,options_int_use,options_float_use)
