# Python 2.7 
# This draft: June 8th, 2014
#
# ----------------------------- Python Configuration --------------------------
from __future__ import division
import numpy as np
import time
#import pylab
#import matplotlib.pyplot as plt
#import matplotlib.mlab as mlab
#import os
#from statsmodels.tsa.filters import hpfilter as hp
import scipy.stats
from scipy import random
from scipy.stats.stats import pearsonr
from scipy.stats import scoreatpercentile as perc
from scipy.optimize import fmin, anneal, brute

# ----------------------------- Preliminaries --------------------------
### Load up the data
T = 72
lambda_ule = np.genfromtxt('lambda_ule_output_2007.csv',delimiter=',').flatten()[:T]
lambda_use = np.genfromtxt('lambda_use_output_2007.csv',delimiter=',').flatten()[:T]
lambda_ne = np.genfromtxt('lambda_ne_output_2007.csv',delimiter=',').flatten()[:T]
u_short = np.genfromtxt('U_shortterm_2007.csv',delimiter=',').flatten()[:T]
u_long = np.genfromtxt('U_longterm_2007.csv',delimiter=',').flatten()[:T]
vac = np.genfromtxt('vac_2007.csv',delimiter=',').flatten()[:T]
N = np.genfromtxt('N_tot_2007.csv',delimiter=',')[:T].flatten()

# -------------------------- Function -----------------------------------
def residual_original(parms):
    # Unpack
    s, m_0, alpha = parms[0], parms[1], parms[2]
    # Make equations
    Eq_1 = lambda_ue - m_0 * A_bar * (vac / (u_short + u_long + s * N))**(1 - alpha)
    Eq_2 = lambda_ne - m_0 * s * (vac / (u_short + u_long + s * N))**(1 - alpha)
    return float((Eq_1**2 + Eq_2**2).sum())

def residual_new(parms):
    # Unpack
    s, m_0, alpha, gamma_1 = parms[0], parms[1], parms[2], parms[3]
    # Make equations
    Eq_1 = lambda_use - m_0 * (vac / (u_short + gamma_1 * u_long + s * N))**(1 - alpha)
    Eq_2 = lambda_ule - m_0 * gamma_1 * (vac / (u_short + gamma_1 * u_long + s * N))**(1 - alpha)
    Eq_3 = lambda_ne - m_0 * s * (vac / (u_short + gamma_1 * u_long + s * N))**(1 - alpha)
    return float((Eq_1**2 + Eq_2**2 + Eq_3**2).sum())

x0 = fmin(residual_new,np.array([.2,.5,.3,.6]))

# -------------------------- Bootstrap -----------------------------------
def residual_boot(parms,lambda_ule_boot,lambda_use_boot,lambda_ne_boot,u_short_boot,u_long_boot,vac_boot,N_boot):
    # Unpack
    s, m_0, alpha, gamma_1 = parms[0], parms[1], parms[2], parms[3]
    # Make equations
    Eq_1 = lambda_use_boot - m_0 * (vac_boot / (u_short_boot + gamma_1 * u_long_boot + s * N_boot))**(1 - alpha)
    Eq_2 = lambda_ule_boot - m_0 * gamma_1 * (vac_boot / (u_short_boot + gamma_1 * u_long_boot + s * N_boot))**(1 - alpha)
    Eq_3 = lambda_ne_boot - m_0 * s * (vac_boot / (u_short_boot + gamma_1 * u_long_boot + s * N_boot))**(1 - alpha)
    return float((Eq_1**2 + Eq_2**2 + Eq_3**2).sum())

N_bootstrap = 2500

# Preallocate

xBootstrap = np.zeros((N_bootstrap,4))

# Bootstrap
start = time.time()
for i in range(N_bootstrap):

    lambda_ule_boot, lambda_use_boot, lambda_ne_boot, u_short_boot, u_long_boot, vac_boot, N_boot = np.zeros((T,1)), np.zeros((T,1)), np.zeros((T,1)), np.zeros((T,1)), np.zeros((T,1)), np.zeros((T,1)), np.zeros((T,1))
    
    for t in range(T):
    
        t_hat = random.randint(0,T-1)
        lambda_ule_boot[t] = lambda_ule[t_hat]
        lambda_use_boot[t] = lambda_use[t_hat]
        lambda_ne_boot[t] = lambda_ne[t_hat]
        u_short_boot[t] = u_short[t_hat]
        u_long_boot[t] = u_long[t_hat]
        vac_boot[t] = vac[t_hat]
        N_boot[t] = N[t_hat]
        
    xBootstrap[i,:] = fmin(residual_boot,x0,args=(lambda_ule_boot,lambda_use_boot,lambda_ne_boot,u_short_boot,u_long_boot,vac_boot,N_boot),disp=0)
end = time.time()


print 'Time ellaspsed', (end - start)

print 'Point estimate', x0

print 'Mean of s', np.mean(xBootstrap[:,0])
print 'Mean of m_0', np.mean(xBootstrap[:,1])
print 'Mean of alpha', np.mean(xBootstrap[:,2])
print 'Mean of gamma_1', np.mean(xBootstrap[:,3])

print 'Standard deviation of s', np.std(xBootstrap[:,0])
print 'Standard deviation of m_0', np.std(xBootstrap[:,1])
print 'Standard deviation of alpha', np.std(xBootstrap[:,2])
print 'Standard deviation of gamma_1', np.std(xBootstrap[:,3])
