import numpy as np
from scipy import stats

import statsmodels.api as sm
from statsmodels.miscmodels import TLinearModel
from statsmodels.tools.numdiff import approx_hess

# Example:
# np.random.seed(98765678)
nobs = 50
nvars = 6
df = 3
rvs = np.random.randn(nobs, nvars - 1)
data_exog = sm.add_constant(rvs, prepend=False)
xbeta = 0.9 + 0.1 * rvs.sum(1)
data_endog = xbeta + 0.1 * np.random.standard_t(df, size=nobs)
print("variance of endog:", data_endog.var())
print("true parameters:", [0.1] * nvars + [0.9])

res_ols = sm.OLS(data_endog, data_exog).fit()
print("\nResults with ols")
print("----------------")
print(res_ols.scale)
print(np.sqrt(res_ols.scale))
print(res_ols.params)
print(res_ols.bse)
kurt = stats.kurtosis(res_ols.resid)
df_fromkurt = 6.0 / kurt + 4
print("df_fromkurt from ols residuals", df_fromkurt)
print(stats.t.stats(df_fromkurt, moments="mvsk"))
print(stats.t.stats(df, moments="mvsk"))

modp = TLinearModel(data_endog, data_exog)
start_value = 0.1 * np.ones(data_exog.shape[1] + 2)
# start_value = np.zeros(data_exog.shape[1]+2)
# start_value[:nvars] = sm.OLS(data_endog, data_exog).fit().params
start_value[:nvars] = res_ols.params
start_value[-2] = df_fromkurt  # 10
start_value[-1] = np.sqrt(res_ols.scale)  # 0.5
modp.start_params = start_value

# adding fixed parameters

fixdf = np.nan * np.zeros(modp.start_params.shape)
fixdf[-2] = 5

fixone = 0
if fixone:
    modp.fixed_params = fixdf
    modp.fixed_paramsmask = np.isnan(fixdf)
    modp.start_params = modp.start_params[modp.fixed_paramsmask]
else:
    modp.fixed_params = None
    modp.fixed_paramsmask = None


print("\nResults with TLinearModel")
print("-------------------------")
resp = modp.fit(
    start_params=modp.start_params, disp=1, method="nm", maxfun=10000, maxiter=5000
)  # 'newton')
# resp = modp.fit(start_params = modp.start_params, disp=1, method='newton')

print("using Nelder-Mead")
print(resp.params)
print(resp.bse)
resp2 = modp.fit(start_params=resp.params, method="Newton")
print("using Newton")
print(resp2.params)
print(resp2.bse)


hb = -approx_hess(modp.start_params, modp.loglike, epsilon=-1e-4)
tmp = modp.loglike(modp.start_params)
print(tmp.shape)
print("eigenvalues of numerical Hessian")
print(np.linalg.eigh(np.linalg.inv(hb))[0])

# store_params is only available in original test script
# pp=np.array(store_params)
# print pp.min(0)
# print pp.max(0)
