import numpy as np
#import scikits.statsmodels.linear_model.regression as smreg

from scipy import special

import scikits.statsmodels.api as sm
from scikits.statsmodels.datasets.macrodata import data
dta = data.load()
gdp = np.log(dta.data['realgdp'])

from numpy import polynomial
from scipy import special


maxorder = 20
polybase = special.chebyt
polybase = special.legendre

t = np.linspace(-1,1,len(gdp))

exog = np.column_stack([polybase(i)(t) for i in range(maxorder)])

fitted = [sm.OLS(gdp, exog[:, :maxr]).fit().fittedvalues for maxr in
          range(2,maxorder)]

print (np.corrcoef(exog[:,1:6], rowvar=0)*10000).astype(int)

import matplotlib.pyplot as plt

plt.figure()
plt.plot(gdp, 'o')
for i in range(maxorder-2):
    plt.plot(fitted[i])

plt.figure()
#plt.plot(gdp, 'o')
for i in range(maxorder-4, maxorder-2):
    #plt.figure()
    plt.plot(gdp - fitted[i])
    plt.title(str(i+2))

plt.figure()
plt.plot(gdp, '.')
plt.plot(fitted[-1], lw=2, color='r')
plt.plot(fitted[0], lw=2, color='g')
plt.title('GDP and Polynomial Trend')

plt.figure()
plt.plot(gdp - fitted[-1], lw=2, color='r')
plt.plot(gdp - fitted[0], lw=2, color='g')
plt.title('Residual GDP minus Polynomial Trend (green: linear, red: legendre(20))')


#orthonormalize an exog using QR

ex2 = t[:,None]**np.arange(6)  #np.vander has columns reversed
q2,r2 = np.linalg.qr(ex2, mode='full')
np.max(np.abs(np.dot(q2.T, q2)-np.eye(6)))
plt.figure()
plt.plot(q2, lw=2)


plt.show()