# -*- coding: utf-8 -*-
"""
Created on Fri Nov 04 10:51:39 2011
@author: josef
"""
import numpy as np
from statsmodels.sandbox.nonparametric import smoothers
from statsmodels.regression.linear_model import OLS, WLS
#DGP: simple polynomial
order = 3
sigma_noise = 0.5
nobs = 100
lb, ub = -1, 2
x = np.linspace(lb, ub, nobs)
x = np.sin(x)
exog = x[:,None]**np.arange(order+1)
y_true = exog.sum(1)
y = y_true + sigma_noise * np.random.randn(nobs)
#xind = np.argsort(x)
pmod = smoothers.PolySmoother(2, x)
pmod.fit(y) #no return
y_pred = pmod.predict(x)
error = y - y_pred
mse = (error*error).mean()
print(mse)
res_ols = OLS(y, exog[:,:3]).fit()
print(np.squeeze(pmod.coef) - res_ols.params)
weights = np.ones(nobs)
weights[:nobs//3] = 0.1
weights[-nobs//5:] = 2
pmodw = smoothers.PolySmoother(2, x)
pmodw.fit(y, weights=weights) #no return
y_predw = pmodw.predict(x)
error = y - y_predw
mse = (error*error).mean()
print(mse)
res_wls = WLS(y, exog[:,:3], weights=weights).fit()
print(np.squeeze(pmodw.coef) - res_wls.params)
doplot = 1
if doplot:
import matplotlib.pyplot as plt
plt.plot(y, '.')
plt.plot(y_true, 'b-', label='true')
plt.plot(y_pred, '-', label='poly')
plt.plot(y_predw, '-', label='poly -w')
plt.legend(loc='upper left')
plt.close()
#plt.show()