Notebook

Linear Regression¶

Linear regression is a technique that measures the relationship between two variables. If we have an independent variable $X$, and a dependent outcome variable $Y$, linear regression allows us to determine which linear model $Y = \alpha + \beta X$ best explains the data.

Python's statsmodels library has a built-in linear fit function. Note that this will give a line of best fit; whether or not the relationship it shows is significant is for you to determine. The output will also have some statistics about the model, such as R-squared and the F value, which may help you quantify how good the fit actually is.

In [50]:
# Import libraries
import numpy as np
from statsmodels import regression
import statsmodels.api as sm
import matplotlib.pyplot as plt
import math
In [65]:
regression.linear_model.OLS?

First we'll define a function that performs linear regression and plots the results.

In [2]:
def linreg(X,Y):
    # Running the linear regression
    # ask the model for alpha too
    X = sm.add_constant(X)
    model = regression.linear_model.OLS(Y, X).fit()
    # alpha
    a = model.params[0] 
    # beta
    b = model.params[1]
    X = X[:, 1]

    # Return summary of the regression and plot results
    X2 = np.linspace(X.min(), X.max(), 100)
    Y_hat = X2 * b + a
    plt.scatter(X, Y, alpha=0.3) # Plot the raw data
    plt.plot(X2, Y_hat, 'r', alpha=0.9);  # Add the regression line, colored in red
    plt.xlabel('X Value')
    plt.ylabel('Y Value')
    return model.summary()
In [48]:
start = '2017-06-01'
end = '2017-07-01'
asset1 = get_pricing('SPY', fields='price', start_date=start, end_date=end, frequency='minute')
ratio1 = 1
asset2 = get_pricing('SPXU', fields='price', start_date=start, end_date=end, frequency='minute')
ratio2 = 16
assetPair = asset1*ratio1-asset2*ratio2
benchmark = get_pricing('SNP', fields='price', start_date=start, end_date=end, frequency='minute')
In [56]:
assetPair.describe()
Out[56]:
count    8580.000000
mean       -5.027272
std         3.146530
min       -15.989000
25%        -7.098500
50%        -4.748000
75%        -2.915000
max         2.835000
dtype: float64
In [57]:
benchmark.describe()
Out[57]:
count    8580.000000
mean       81.108227
std         1.622526
min        77.740000
25%        79.610000
50%        81.510000
75%        82.630000
max        83.600000
Name: Equity(22169 [SNP]), dtype: float64
In [61]:
# We have to take the percent changes to get to returns
# Get rid of the first (0th) element because it is NAN
pct_change_asset = assetPair.pct_change()[1:]
pct_change_benchmark = benchmark.pct_change()[1:]
In [62]:
pct_change_asset.describe()
Out[62]:
count    8579.000000
mean        0.004704
std         0.384204
min        -7.625000
25%        -0.018783
50%         0.000000
75%         0.019686
max        30.333333
dtype: float64
In [63]:
pct_change_benchmark.describe()
Out[63]:
count    8579.000000
mean       -0.000005
std         0.000532
min        -0.015716
25%         0.000000
50%         0.000000
75%         0.000000
max         0.010351
Name: Equity(22169 [SNP]), dtype: float64
In [68]:
# Y axis: asset change
# X axis: benchmark change
linreg(pct_change_asset.values, pct_change_benchmark.values)
Out[68]:
<caption>OLS Regression Results</caption>
Dep. Variable: y R-squared: 0.001
Model: OLS Adj. R-squared: 0.001
Method: Least Squares F-statistic: 5.862
Date: Mon, 03 Jul 2017 Prob (F-statistic): 0.0155
Time: 04:49:46 Log-Likelihood: 52512.
No. Observations: 8579 AIC: -1.050e+05
Df Residuals: 8577 BIC: -1.050e+05
Df Model: 1
Covariance Type: nonrobust
coef std err t P>|t| [95.0% Conf. Int.]
const -5.194e-06 5.74e-06 -0.905 0.365 -1.64e-05 6.05e-06
x1 -3.616e-05 1.49e-05 -2.421 0.015 -6.54e-05 -6.88e-06
Omnibus: 10180.265 Durbin-Watson: 2.188
Prob(Omnibus): 0.000 Jarque-Bera (JB): 19102941.885
Skew: -5.274 Prob(JB): 0.00
Kurtosis: 233.932 Cond. No. 2.60

Each point on the above graph represents a day, with the x-coordinate being the return of SPY, and the y-coordinate being the return of TSLA. As we can see, the line of best fit tells us that for every 1% increased return we see from the SPY, we should see an extra 1.92% from TSLA. This is expressed by the parameter $\beta$, which is 1.9271 as estimated. Of course, for decresed return we will also see about double the loss in TSLA, so we haven't gained anything, we are just more volatile.

Linear Regression vs. Correlation¶

  • Linear regression gives us a specific linear model, but is limited to cases of linear dependence.
  • Correlation is general to linear and non-linear dependencies, but doesn't give us an actual model.
  • Both are measures of covariance.
  • Linear regression can give us relationship between Y and many independent variables by making X multidimensional.

Knowing Parameters vs. Estimates¶

It is very important to keep in mind that all $\alpha$ and $\beta$ parameters estimated by linear regression are just that - estimates. You can never know the underlying true parameters unless you know the physical process producing the data. The paremeters you estimate today may not be the same as the same analysis done including tomorrow's data, and the underlying true parameters may be moving. As such it is very important when doing actual analysis to pay attention to the standard error of the parameter estimates. More material on the standard error will be presented in a later lecture. One way to get a sense of how stable your paremeter estimates are is to estimates them using a rolling window of data and see how much variance there is in the estimates.

Evaluating and reporting results¶

The regression model relies on several assumptions:

  • The independent variable is not random.
  • The variance of the error term is constant across observations. This is important for evaluating the goodness of the fit.
  • The errors are not autocorrelated. The Durbin-Watson statistic detects this; if it is close to 2, there is no autocorrelation.
  • The errors are normally distributed. If this does not hold, we cannot use some of the statistics, such as the F-test.

If we confirm that the necessary assumptions of the regression model are satisfied, we can safely use the statistics reported to analyze the fit. For example, the $R^2$ value tells us the fraction of the total variation of $Y$ that is explained by the model.

When making a prediction based on the model, it's useful to report not only a single value but a confidence interval. The linear regression reports 95% confidence intervals for the regression parameters, and we can visualize what this means using the seaborn library, which plots the regression line and highlights the 95% (by default) confidence interval for the regression line:

In [69]:
import seaborn

seaborn.regplot(pct_change_asset.values, pct_change_benchmark.values);

Mathematical Background¶

This is a very brief overview of linear regression. For more, please see: https://en.wikipedia.org/wiki/Linear_regression

Ordinary Least Squares¶

Regression works by optimizing the placement of the line of best fit (or plane in higher dimensions). It does so by defining how bad the fit is using an objective function. In ordinary least squares regression (OLS), what we use here, the objective function is:

$$\sum_{i=1}^n (Y_i - a - bX_i)^2$$

We use $a$ and $b$ to represent the potential candidates for $\alpha$ and $\beta$. What this objective function means is that for each point on the line of best fit we compare it with the real point and take the square of the difference. This function will decrease as we get better parameter estimates. Regression is a simple case of numerical optimization that has a closed form solution and does not need any optimizer. We just find the results that minimize the objective function.

We will denote the eventual model that results from minimizing our objective function as:

$$ \hat{Y} = \hat{\alpha} + \hat{\beta}X $$

With $\hat{\alpha}$ and $\hat{\beta}$ being the chosen estimates for the parameters that we use for prediction and $\hat{Y}$ being the predicted values of $Y$ given the estimates.

Standard Error¶

We can also find the standard error of estimate, which measures the standard deviation of the error term $\epsilon$, by getting the scale parameter of the model returned by the regression and taking its square root. The formula for standard error of estimate is $$ s = \left( \frac{\sum_{i=1}^n \epsilon_i^2}{n-2} \right)^{1/2} $$

If $\hat{\alpha}$ and $\hat{\beta}$ were the true parameters ($\hat{\alpha} = \alpha$ and $\hat{\beta} = \beta$), we could represent the error for a particular predicted value of $Y$ as $s^2$ for all values of $X_i$. We could simply square the difference $(Y - \hat{Y})$ to get the variance because $\hat{Y}$ incorporates no error in the paremeter estimates themselves. Because $\hat{\alpha}$ and $\hat{\beta}$ are merely estimates in our construction of the model of $Y$, any predicted values , $\hat{Y}$, will have their own standard error based on the distribution of the $X$ terms that we plug into the model. This forecast error is represented by the following:

$$ s_f^2 = s^2 \left( 1 + \frac{1}{n} + \frac{(X - \mu_X)^2}{(n-1)\sigma_X^2} \right) $$

where $\mu_X$ is the mean of our observations of $X$ and $\sigma_X$ is the standard deviation of $X$. This adjustment to $s^2$ incorporates the uncertainty in our parameter estimates. Then the 95% confidence interval for the prediction is $\hat{Y} \pm t_cs_f$, where $t_c$ is the critical value of the t-statistic for $n$ samples and a desired 95% confidence.