<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Rolling Regression
Author
Step1:  Lets load the prices of GFI and GLD.
Step2:  Plotting the prices over time suggests a strong correlation. However, the correlation seems to change over time.
Step3:  A naive approach would be to estimate a linear model and ignore the time domain.
Step4:  The posterior predictive plot shows how bad the fit is.
Step5:  Rolling regression
Next, we will build an improved model that will allow for changes in the regression coefficients over time. Specifically, we will assume that intercept and slope follow a random-walk through time. That idea is similar to the stochastic volatility model.
$$ \alpha_t \sim \mathcal{N}(\alpha_{t-1}, \sigma_\alpha^2) $$
$$ \beta_t \sim \mathcal{N}(\beta_{t-1}, \sigma_\beta^2) $$
First, lets define the hyper-priors for $\sigma_\alpha^2$ and $\sigma_\beta^2$. This parameter can be interpreted as the volatility in the regression coefficients.
Step6:  Next, we define the regression parameters that are not a single random variable but rather a random vector with the above stated dependence structure. So as not to fit a coefficient to a single data point, we will chunk the data into bins of 50 and apply the same coefficients to all data points in a single bin.
Step7:  Perform the regression given coefficients and data and link to the data via the likelihood.
Step8:  Inference. Despite this being quite a complex model, NUTS handles it wells.
Step9:  Analysis of results
$\alpha$, the intercept, does not seem to change over time.
Step10:  However, the slope does.
Step11:  The posterior predictive plot shows that we capture the change in regression over time much better. Note that we should have used returns instead of prices. The model would still work the same, but the visualisations would not be quite as clear.
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
import pandas as pd
import numpy as np
import pymc3 as pm
import matplotlib.pyplot as plt
Explanation: Rolling Regression
Author: Thomas Wiecki
Pairs trading is a famous technique in algorithmic trading that plays two stocks against each other.
For this to work, stocks must be correlated (cointegrated).
One common example is the price of gold (GLD) and the price of gold mining operations (GFI).
End of explanation
# from pandas_datareader import data
# prices = data.GoogleDailyReader(symbols=['GLD', 'GFI'], end='2014-8-1').read().loc['Open', :, :]
prices = pd.read_csv(pm.get_data('stock_prices.csv'))
prices['Date'] = pd.DatetimeIndex(prices['Date'])
prices = prices.set_index('Date')
prices.head()
finite_idx = (np.isfinite(prices.GLD.values)) & (np.isfinite(prices.GFI.values))
prices = prices.iloc[finite_idx]
Explanation: Lets load the prices of GFI and GLD.
End of explanation
fig = plt.figure(figsize=(9, 6))
ax = fig.add_subplot(111, xlabel='Price GFI in \$', ylabel='Price GLD in \$')
colors = np.linspace(0.1, 1, len(prices))
mymap = plt.get_cmap("winter")
sc = ax.scatter(prices.GFI, prices.GLD, c=colors, cmap=mymap, lw=0)
cb = plt.colorbar(sc)
cb.ax.set_yticklabels([str(p.date()) for p in prices[::len(prices)//10].index]);
Explanation: Plotting the prices over time suggests a strong correlation. However, the correlation seems to change over time.
End of explanation
with pm.Model() as model_reg:
    pm.glm.GLM.from_formula('GLD ~ GFI', prices)
    trace_reg = pm.sample(2000)
Explanation: A naive approach would be to estimate a linear model and ignore the time domain.
End of explanation
fig = plt.figure(figsize=(9, 6))
ax = fig.add_subplot(111, xlabel='Price GFI in \$', ylabel='Price GLD in \$', 
            title='Posterior predictive regression lines')
sc = ax.scatter(prices.GFI, prices.GLD, c=colors, cmap=mymap, lw=0)
pm.plot_posterior_predictive_glm(trace_reg[100:], samples=100, 
                              label='posterior predictive regression lines',
                              lm=lambda x, sample: sample['Intercept'] + sample['GFI'] * x,
                              eval=np.linspace(prices.GFI.min(), prices.GFI.max(), 100))
cb = plt.colorbar(sc)
cb.ax.set_yticklabels([str(p.date()) for p in prices[::len(prices)//10].index]);
ax.legend(loc=0);
Explanation: The posterior predictive plot shows how bad the fit is.
End of explanation
model_randomwalk = pm.Model()
with model_randomwalk:
    # std of random walk, best sampled in log space.
    sigma_alpha = pm.Exponential('sigma_alpha', 1./.02, testval = .1)
    sigma_beta = pm.Exponential('sigma_beta', 1./.02, testval = .1)
Explanation: Rolling regression
Next, we will build an improved model that will allow for changes in the regression coefficients over time. Specifically, we will assume that intercept and slope follow a random-walk through time. That idea is similar to the stochastic volatility model.
$$ \alpha_t \sim \mathcal{N}(\alpha_{t-1}, \sigma_\alpha^2) $$
$$ \beta_t \sim \mathcal{N}(\beta_{t-1}, \sigma_\beta^2) $$
First, lets define the hyper-priors for $\sigma_\alpha^2$ and $\sigma_\beta^2$. This parameter can be interpreted as the volatility in the regression coefficients.
End of explanation
import theano.tensor as tt
# To make the model simpler, we will apply the same coefficient for 50 data points at a time
subsample_n = 50
lendata = len(prices)
ncoef = lendata // subsample_n
idx = range(ncoef * subsample_n)
with model_randomwalk:
    alpha = pm.GaussianRandomWalk('alpha', sigma_alpha**-2, 
                                  shape=ncoef)
    beta = pm.GaussianRandomWalk('beta', sigma_beta**-2, 
                                 shape=ncoef)
    
    # Make coefficients have the same length as prices
    alpha_r = tt.repeat(alpha, subsample_n)
    beta_r = tt.repeat(beta, subsample_n)    
Explanation: Next, we define the regression parameters that are not a single random variable but rather a random vector with the above stated dependence structure. So as not to fit a coefficient to a single data point, we will chunk the data into bins of 50 and apply the same coefficients to all data points in a single bin.
End of explanation
with model_randomwalk:
    # Define regression
    regression = alpha_r + beta_r * prices.GFI.values[idx]
    
    # Assume prices are Normally distributed, the mean comes from the regression.
    sd = pm.Uniform('sd', 0, 20)
    likelihood = pm.Normal('y', 
                           mu=regression, 
                           sd=sd, 
                           observed=prices.GLD.values[idx])
Explanation: Perform the regression given coefficients and data and link to the data via the likelihood.
End of explanation
with model_randomwalk:
    trace_rw = pm.sample(2000, njobs=2)
Explanation: Inference. Despite this being quite a complex model, NUTS handles it wells.
End of explanation
fig = plt.figure(figsize=(8, 6))
ax = plt.subplot(111, xlabel='time', ylabel='alpha', title='Change of alpha over time.')
ax.plot(trace_rw[-1000:]['alpha'].T, 'r', alpha=.05);
ax.set_xticklabels([str(p.date()) for p in prices[::len(prices)//5].index]);
Explanation: Analysis of results
$\alpha$, the intercept, does not seem to change over time.
End of explanation
fig = plt.figure(figsize=(8, 6))
ax = fig.add_subplot(111, xlabel='time', ylabel='beta', title='Change of beta over time')
ax.plot(trace_rw[-1000:]['beta'].T, 'b', alpha=.05);
ax.set_xticklabels([str(p.date()) for p in prices[::len(prices)//5].index]);
Explanation: However, the slope does.
End of explanation
fig = plt.figure(figsize=(8, 6))
ax = fig.add_subplot(111, xlabel='Price GFI in \$', ylabel='Price GLD in \$', 
            title='Posterior predictive regression lines')
colors = np.linspace(0.1, 1, len(prices))
colors_sc = np.linspace(0.1, 1, len(trace_rw[-500::10]['alpha'].T))
mymap = plt.get_cmap('winter')
mymap_sc = plt.get_cmap('winter')
xi = np.linspace(prices.GFI.min(), prices.GFI.max(), 50)
for i, (alpha, beta) in enumerate(zip(trace_rw[-500::10]['alpha'].T, trace_rw[-500::10]['beta'].T)):
    for a, b in zip(alpha, beta):
        ax.plot(xi, a + b*xi, alpha=.05, lw=1, c=mymap_sc(colors_sc[i]))
        
sc = ax.scatter(prices.GFI, prices.GLD, label='data', cmap=mymap, c=colors)
cb = plt.colorbar(sc)
cb.ax.set_yticklabels([str(p.date()) for p in prices[::len(prices)//10].index]);
Explanation: The posterior predictive plot shows that we capture the change in regression over time much better. Note that we should have used returns instead of prices. The model would still work the same, but the visualisations would not be quite as clear.
End of explanation
<END_TASK>