# -*- coding: utf-8 -*-
"""
Created on Fri Feb 12 01:01:50 2010

Author: josef-pktd

latest result
-------------
all are very close
garch0 has different parameterization of constant
ordering of parameters is different


seed 2780185
h.shape (2000,)
Optimization terminated successfully.
         Current function value: 2093.813397
         Iterations: 387
         Function evaluations: 676
ggres.params [-0.6146253   0.1914537   0.01039355  0.78802188]
Optimization terminated successfully.
         Current function value: 2093.972953
         Iterations: 201
         Function evaluations: 372
ggres0.params [-0.61537527  0.19635128  4.00706058]
Warning: Desired error not necessarily achieveddue to precision loss
         Current function value: 2093.972953
         Iterations: 51
         Function evaluations: 551
         Gradient evaluations: 110
ggres0.params [-0.61537855  0.19635265  4.00694669]
Optimization terminated successfully.
         Current function value: 2093.751420
         Iterations: 103
         Function evaluations: 187
[ 0.78671519  0.19692222  0.61457171]
-2093.75141963

Final Estimate:
 LLH:  2093.750    norm LLH:  2.093750 
    omega    alpha1     beta1 
0.7867438 0.1970437 0.6145467 

long run variance comparison
^^^^^^^^^^^^^^^^^^^^^^^^^^^^
R
>>> 0.7867438/(1- 0.1970437- 0.6145467)
4.1757097302897526
Garch (gjr) asymetric, longrun var ?
>>> 1/(1-0.6146253 - 0.1914537 - 0.01039355) * 0.78802188
4.2937548579245242
>>> 1/(1-0.6146253 - 0.1914537 + 0.01039355) * 0.78802188
3.8569053452140345
Garch0
>>> (1-0.61537855 - 0.19635265) * 4.00694669
0.7543830449902722
>>> errgjr4.var() #for different random seed
4.0924199964716106

todo: add code and verify, check for longer lagpolys

"""


import numpy as np
from numpy.testing import assert_almost_equal

import matplotlib.pyplot as plt
import numdifftools as ndt

import scikits.statsmodels.api as sm
from scikits.statsmodels.sandbox import tsa
from scikits.statsmodels.sandbox.tsa.garch import *  # local import


nobs = 1000
examples = ['garch', 'rpyfit']
if 'garch' in examples:
    err,h = generate_kindofgarch(nobs, [1.0, -0.95], [1.0,  0.1], mu=0.5)
    plt.figure()
    plt.subplot(211)
    plt.plot(err)
    plt.subplot(212)
    plt.plot(h)
    #plt.show()
    
    seed = 3842774 #91234  #8837708
    seed = np.random.randint(9999999)
    print 'seed', seed
    np.random.seed(seed)
    ar1 = -0.9
    err,h = generate_garch(nobs, [1.0, ar1], [1.0,  0.50], mu=0.0,scale=0.1)
#    plt.figure()
#    plt.subplot(211)
#    plt.plot(err)
#    plt.subplot(212)
#    plt.plot(h)
#    plt.figure()
#    plt.subplot(211)
#    plt.plot(err[-400:])
#    plt.subplot(212)
#    plt.plot(h[-400:])
    #plt.show()
    garchplot(err, h)
    garchplot(err[-400:], h[-400:])
 
 
    np.random.seed(seed)
    errgjr,hgjr, etax = generate_gjrgarch(nobs, [1.0, ar1], 
                                [[1,0],[0.5,0]], mu=0.0,scale=0.1)    
    garchplot(errgjr[:nobs], hgjr[:nobs], 'GJR-GARCH(1,1) Simulation - symmetric')
    garchplot(errgjr[-400:nobs], hgjr[-400:nobs], 'GJR-GARCH(1,1) Simulation - symmetric')
    
    np.random.seed(seed)
    errgjr2,hgjr2, etax = generate_gjrgarch(nobs, [1.0, ar1], 
                                [[1,0],[0.1,0.9]], mu=0.0,scale=0.1)    
    garchplot(errgjr2[:nobs], hgjr2[:nobs], 'GJR-GARCH(1,1) Simulation')
    garchplot(errgjr2[-400:nobs], hgjr2[-400:nobs], 'GJR-GARCH(1,1) Simulation')
    
    np.random.seed(seed)
    errgjr3,hgjr3, etax3 = generate_gjrgarch(nobs, [1.0, ar1], 
                        [[1,0],[0.1,0.9],[0.1,0.9],[0.1,0.9]], mu=0.0,scale=0.1)    
    garchplot(errgjr3[:nobs], hgjr3[:nobs], 'GJR-GARCH(1,3) Simulation')
    garchplot(errgjr3[-400:nobs], hgjr3[-400:nobs], 'GJR-GARCH(1,3) Simulation')

    np.random.seed(seed)
    errgjr4,hgjr4, etax4 = generate_gjrgarch(nobs, [1.0, ar1], 
                        [[1., 1,0],[0, 0.1,0.9],[0, 0.1,0.9],[0, 0.1,0.9]], 
                        mu=0.0,scale=0.1)    
    garchplot(errgjr4[:nobs], hgjr4[:nobs], 'GJR-GARCH(1,3) Simulation')
    garchplot(errgjr4[-400:nobs], hgjr4[-400:nobs], 'GJR-GARCH(1,3) Simulation')
    
    varinno = np.zeros(100)
    varinno[0] = 1.
    errgjr5,hgjr5, etax5 = generate_gjrgarch(100, [1.0, -0.], 
                        [[1., 1,0],[0, 0.1,0.8],[0, 0.05,0.7],[0, 0.01,0.6]], 
                        mu=0.0,scale=0.1, varinnovation=varinno)
    garchplot(errgjr5[:20], hgjr5[:20], 'GJR-GARCH(1,3) Simulation')
    #garchplot(errgjr4[-400:nobs], hgjr4[-400:nobs], 'GJR-GARCH(1,3) Simulation')
    

#plt.show()
seed = np.random.randint(9999999)  # 9188410
print 'seed', seed
    
x = np.arange(20).reshape(10,2)
x3 = np.column_stack((np.ones((x.shape[0],1)),x))
y, inp = miso_lfilter([1., 0],np.array([[-2.0,3,1],[0.0,0.0,0]]),x3)

nobs = 1000
warmup = 1000
np.random.seed(seed)
ar = [1.0, -0.7]#7, -0.16, -0.1]
#ma = [[1., 1, 0],[0, 0.6,0.1],[0, 0.1,0.1],[0, 0.1,0.1]]
ma = [[1., 0, 0],[0, 0.8,0.0]] #,[0, 0.9,0.0]]
#    errgjr4,hgjr4, etax4 = generate_gjrgarch(warmup+nobs, [1.0, -0.99], 
#                        [[1., 1, 0],[0, 0.6,0.1],[0, 0.1,0.1],[0, 0.1,0.1]], 
#                        mu=0.2, scale=0.25)   

errgjr4,hgjr4, etax4 = generate_gjrgarch(warmup+nobs, ar, ma, 
                     mu=0.4, scale=1.01)
errgjr4,hgjr4, etax4 = errgjr4[warmup:], hgjr4[warmup:], etax4[warmup:]
garchplot(errgjr4[:nobs], hgjr4[:nobs], 'GJR-GARCH(1,3) Simulation - DGP')
ggmod = Garch(errgjr4-errgjr4.mean())#hgjr4[:nobs])#-hgjr4.mean()) #errgjr4)
ggmod.nar = 1
ggmod.nma = 1
ggmod._start_params = np.array([-0.6, 0.1, 0.2, 0.0])
ggres = ggmod.fit(start_params=np.array([-0.6, 0.1, 0.2, 0.0]), maxiter=1000)
print 'ggres.params', ggres.params
garchplot(ggmod.errorsest, ggmod.h, title='Garch estimated')

ggmod0 = Garch0(errgjr4-errgjr4.mean())#hgjr4[:nobs])#-hgjr4.mean()) #errgjr4)
ggmod0.nar = 1
ggmod.nma = 1
start_params = np.array([-0.6, 0.2, 0.1])
ggmod0._start_params = start_params #np.array([-0.6, 0.1, 0.2, 0.0])
ggres0 = ggmod0.fit(start_params=start_params, maxiter=2000)
print 'ggres0.params', ggres0.params

ggmod0 = Garch0(errgjr4-errgjr4.mean())#hgjr4[:nobs])#-hgjr4.mean()) #errgjr4)
ggmod0.nar = 1
ggmod.nma = 1
start_params = np.array([-0.6, 0.2, 0.1])
ggmod0._start_params = start_params #np.array([-0.6, 0.1, 0.2, 0.0])
ggres0 = ggmod0.fit(start_params=start_params, method='bfgs', maxiter=2000)
print 'ggres0.params', ggres0.params


g11res = optimize.fmin(lambda params: -loglike_GARCH11(params, errgjr4-errgjr4.mean())[0], [0.93, 0.9, 0.2])
print g11res
llf = loglike_GARCH11(g11res, errgjr4-errgjr4.mean())
print llf[0]

if 'rpyfit' in examples:
    from rpy import r
    r.library('fGarch')
    f = r.formula('~garch(1, 1)')
    fit = r.garchFit(f, data = errgjr4-errgjr4.mean(), include_mean=False)

if 'rpysim' in examples:
    from rpy import r
    f = r.formula('~garch(1, 1)')
    #fit = r.garchFit(f, data = errgjr4)
    x = r.garchSim( n = 500)
    print 'R acf', tsa.acf(np.power(x,2))[:15]
    arma3 = Arma(np.power(x,2))
    arma3res = arma3.fit(start_params=[-0.2,0.1,0.5],maxiter=5000)
    print arma3res.params
    arma3b = Arma(np.power(x,2))
    arma3bres = arma3b.fit(start_params=[-0.2,0.1,0.5],maxiter=5000, method='bfgs')
    print arma3bres.params
    
    xr = r.garchSim( n = 100)
    
    x = np.asarray(xr)
    ggmod = Garch(x-x.mean())
    ggmod.nar = 1
    ggmod.nma = 1
    ggmod._start_params = np.array([-0.6, 0.1, 0.2, 0.0])
    ggres = ggmod.fit(start_params=np.array([-0.6, 0.1, 0.2, 0.0]), maxiter=1000)
    print 'ggres.params', ggres.params
    
    g11res = optimize.fmin(lambda params: -loglike_GARCH11(params, x-x.mean())[0], [0.6, 0.6, 0.2])
    print g11res
    llf = loglike_GARCH11(g11res, x-x.mean())
    print llf[0]
    
    garchplot(ggmod.errorsest, ggmod.h, title='Garch estimated')
    fit = r.garchFit(f, data = x-x.mean(), include_mean=False, trace=False)
    print r.summary(fit)

'''based on R default simulation
model = list(omega = 1e-06, alpha = 0.1, beta = 0.8)
nobs = 1000
(with nobs=500, gjrgarch doesn't do well

>>> ggres = ggmod.fit(start_params=np.array([-0.6, 0.1, 0.2, 0.0]), maxiter=1000)
Optimization terminated successfully.
         Current function value: -448.861335
         Iterations: 385
         Function evaluations: 690
>>> print 'ggres.params', ggres.params
ggres.params [ -7.75090330e-01   1.57714749e-01  -9.60223930e-02   8.76021411e-07]
rearranged 
8.76021411e-07 1.57714749e-01(-9.60223930e-02) 7.75090330e-01

>>> print g11res
[  2.97459808e-06   7.83128600e-01   2.41110860e-01]
>>> llf = loglike_GARCH11(g11res, x-x.mean())
>>> print llf[0]
442.603541936

Log Likelihood:
 -448.9376    normalized:  -4.489376 
      omega       alpha1        beta1  
1.01632e-06  1.02802e-01  7.57537e-01 
'''


''' the following is for errgjr4-errgjr4.mean()
ggres.params [-0.54510407  0.22723132  0.06482633  0.82325803]
Final Estimate:
 LLH:  2065.56    norm LLH:  2.06556 
        mu      omega     alpha1      beta1 
0.07229732 0.83069480 0.26313883 0.53986167 

ggres.params [-0.50779163  0.2236606   0.00700036  1.154832  
Final Estimate:
 LLH:  2116.084    norm LLH:  2.116084 
           mu         omega        alpha1         beta1 
-4.759227e-17  1.145404e+00  2.288348e-01  5.085949e-01

run3
DGP
0.4/??    0.8      0.7
gjrgarch:
ggres.params [-0.45196579  0.2569641   0.02201904  1.11942636]
rearranged
const/omega  ma1/alpha1             ar1/beta1
1.11942636 0.2569641(+0.02201904) 0.45196579
g11:
[ 1.10262688  0.26680468  0.45724957]
-2055.73912687
R:
Final Estimate:
 LLH:  2055.738    norm LLH:  2.055738 
           mu         omega        alpha1         beta1 
-1.665226e-17  1.102396e+00  2.668712e-01  4.573224e-01 
fit = r.garchFit(f, data = errgjr4-errgjr4.mean())
rpy.RPy_RException: Error in solve.default(fit$hessian) : 
  Lapack routine dgesv: system is exactly singular

run4
DGP:
mu=0.4, scale=1.01
ma = [[1., 0, 0],[0, 0.8,0.0]], ar = [1.0, -0.7]
maybe something wrong with simulation

gjrgarch
ggres.params [-0.50554663  0.24449867 -0.00521004  1.00796791]
rearranged
1.00796791   0.24449867(-0.00521004)   0.50554663
garch11:
[ 1.01258264  0.24149155  0.50479994]
-2056.3877404
R include_constant=False
Final Estimate:
 LLH:  2056.397    norm LLH:  2.056397 
    omega    alpha1     beta1 
1.0123560 0.2409589 0.5049154 
'''


erro,ho, etaxo = generate_gjrgarch(20, ar, ma, mu=0.04, scale=0.01, 
                  varinnovation = np.ones(20))

if 'sp500' in examples:
    import tabular as tb
    import scikits.timeseries as ts
    
    a = tb.loadSV(r'C:\Josef\work-oth\gspc_table.csv')
    
    s = ts.time_series(a[0]['Close'][::-1],
                dates=ts.date_array(a[0]['Date'][::-1],freq="D"))
    
    sp500 = a[0]['Close'][::-1]
    sp500r = np.diff(np.log(sp500))


#plt.show()