# -*- coding: utf-8 -*-
"""
Created on Thu Apr 29 17:06:47 2021

@author: arnep

Estimation of the extremal index of data generated by an
an AR(1) model with student-t distributed errors and stable coefficient.

Estimators used are the parametric kappa filter, the stable and rolling horizon
estimator based on Cai (2019), and the sliding block estimators of Berghaus &
Bücher (2018).
"""
    
import numpy as np
import matplotlib.pyplot as plt
from methods_simulation import ground_truth_ARt, random_ARt, MSE_miss_handler
from methods_simulation import est_ARt_parametric, est_non_parametr_stable, est_non_parametr_rol_hor
from methods_simulation import est_sliding_blocks_Z, est_sliding_blocks_Y, theta_simulator
from methods_extremal_index import MLE, construct_psi_k2, construct_psi_robust_k2
from methods_extremal_index import non_parametric_day_based

 
## simulation parameters
ph = 0.8
n = 1000
runs = 10
phi= np.ones(n)*ph
horizon_length = 400
dof = 3
func_gt = ground_truth_ARt
func_rand = random_ARt
func_crit = MSE_miss_handler


#object names
name = 'MSE_ARt_stable_phi_' +str(ph)


#estimation parameters
func_non_param = non_parametric_day_based
k = 50
k2 =  int(round(k/n*horizon_length)*1.2)
bn = int(round(n/k)*4)
d_u = 10
func_mle = MLE
func_filter = construct_psi_robust_k2

### est with robust filter
theta_init = [0.8, 0, 0,  3]
func_filter2 = construct_psi_k2
theta_init2 = [ 0.8,  0, 0,  3]

bounds = ((-2,2), #omega
        (-.9999,  .9999),#alpha
        (-.9999,  .9999), #beta
        (2, 10001) #nu
        )


#create simulation object    
sim_MSE_ARt = theta_simulator(name, n, runs, func_gt, func_rand, func_crit, dof =dof, phi=phi )


#estimate parametric model
func_est = est_ARt_parametric
print('Estimating parametric model:')
sim_MSE_ARt.simulate(func_est, name = 'Parametric ARt(1) filter, kappa: 10',
                     func_filter = func_filter2, theta_init =theta_init2,
                     bounds = bounds, func_mle = func_mle)

#estimate non-parametric estimator of Cai (2019) (stable and rolling horizon)
func_est = est_non_parametr_stable
print('Estimating stable non-parametric estimator:')
sim_MSE_ARt.simulate(func_est, name = 'Stable non-parametric estimator, k: ' + str(k),
                     d_u = d_u, k = k, func = func_non_param)
    
func_est = est_non_parametr_rol_hor
print('Estimating rolling horizon non-parametric estimator:')
sim_MSE_ARt.simulate(func_est, name = 'Rolling horizon non-parametric estimator, k: ' + str(k2) +
                     ', horizon length: ' + str(horizon_length),
                     horizon_length= horizon_length, d_u = d_u, k = k2,
                     func = func_non_param)

#estimate sliding blocks estimator of Berghaus & Bucher (2018)
func_est =  est_sliding_blocks_Y
print('Estimating Y-estimator:')
sim_MSE_ARt.simulate(func_est, name = 'Sliding blocks Y-estimator, bn: ' + str(bn), 
                     bn = bn)

func_est = est_sliding_blocks_Z
print('Estimating Z-estimator:')
sim_MSE_ARt.simulate(func_est, name = 'Sliding blocks Z-estimator, bn: ' + str(bn),
                     bn = bn)

#get the summary
sim_MSE_ARt.calc_stats()

#plot all the estimators
for key in sim_MSE_ARt.theta.keys():
    sim_MSE_ARt.plot_theta_gt(key)

#plot ground truth and simulation example
sim_MSE_ARt.plot_gt('ARt(1) model with phi = ' + str(ph) )
sim_MSE_ARt.plot_sim_example('ARt(1) model with phi = ' + str(ph))

#plot autoregressive coefficient
plt.figure(figsize=(10, 6), dpi=200)
plt.grid(color = 'black')
plt.plot(phi)
plt.xticks(np.arange(0, n*1.2, n/5))
plt.yticks(np.arange(0, 1.2, .2))
plt.title('$\phi_t$')
plt.show()

