Linear Regression, Statistical Learning within R

This text contains:

• Simple linear regression
• Multiple linear regression
• Nonlinear terms and Interactions
• Qualitative predictors fix

# Libaries containing example data
library(MASS)
library(ISLR)

Simple linear regression

names(Boston)

##  [1] "crim"    "zn"      "indus"   "chas"    "nox"     "rm"      "age"    
##  [8] "dis"     "rad"     "tax"     "ptratio" "black"   "lstat"   "medv"

## ?Boston
plot(medv ~ lstat, Boston)
fit1 = lm(medv ~ lstat, data = Boston)
fit1

## 
## Call:
## lm(formula = medv ~ lstat, data = Boston)
## 
## Coefficients:
## (Intercept)        lstat  
##       34.55        -0.95

summary(fit1)

## 
## Call:
## lm(formula = medv ~ lstat, data = Boston)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -15.17  -3.99  -1.32   2.03  24.50 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  34.5538     0.5626    61.4   <2e-16 ***
## lstat        -0.9500     0.0387   -24.5   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.22 on 504 degrees of freedom
## Multiple R-squared:  0.544,  Adjusted R-squared:  0.543 
## F-statistic:  602 on 1 and 504 DF,  p-value: <2e-16

abline(fit1, col = "red")

names(fit1)

##  [1] "coefficients"  "residuals"     "effects"       "rank"         
##  [5] "fitted.values" "assign"        "qr"            "df.residual"  
##  [9] "xlevels"       "call"          "terms"         "model"

confint(fit1)

##              2.5 % 97.5 %
## (Intercept) 33.448 35.659
## lstat       -1.026 -0.874

predict(fit1, data.frame(lstat = c(5, 10, 15)), interval = "confidence")

##     fit   lwr   upr
## 1 29.80 29.01 30.60
## 2 25.05 24.47 25.63
## 3 20.30 19.73 20.87

Multiple linear regression

fit2 = lm(medv ~ lstat + age, data = Boston)
summary(fit2)

## 
## Call:
## lm(formula = medv ~ lstat + age, data = Boston)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -15.98  -3.98  -1.28   1.97  23.16 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  33.2228     0.7308   45.46   <2e-16 ***
## lstat        -1.0321     0.0482  -21.42   <2e-16 ***
## age           0.0345     0.0122    2.83   0.0049 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.17 on 503 degrees of freedom
## Multiple R-squared:  0.551,  Adjusted R-squared:  0.549 
## F-statistic:  309 on 2 and 503 DF,  p-value: <2e-16

fit3 = lm(medv ~ ., Boston)
summary(fit3)

## 
## Call:
## lm(formula = medv ~ ., data = Boston)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.594  -2.730  -0.518   1.777  26.199 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.65e+01   5.10e+00    7.14  3.3e-12 ***
## crim        -1.08e-01   3.29e-02   -3.29  0.00109 ** 
## zn           4.64e-02   1.37e-02    3.38  0.00078 ***
## indus        2.06e-02   6.15e-02    0.33  0.73829    
## chas         2.69e+00   8.62e-01    3.12  0.00193 ** 
## nox         -1.78e+01   3.82e+00   -4.65  4.2e-06 ***
## rm           3.81e+00   4.18e-01    9.12  < 2e-16 ***
## age          6.92e-04   1.32e-02    0.05  0.95823    
## dis         -1.48e+00   1.99e-01   -7.40  6.0e-13 ***
## rad          3.06e-01   6.63e-02    4.61  5.1e-06 ***
## tax         -1.23e-02   3.76e-03   -3.28  0.00111 ** 
## ptratio     -9.53e-01   1.31e-01   -7.28  1.3e-12 ***
## black        9.31e-03   2.69e-03    3.47  0.00057 ***
## lstat       -5.25e-01   5.07e-02  -10.35  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.75 on 492 degrees of freedom
## Multiple R-squared:  0.741,  Adjusted R-squared:  0.734 
## F-statistic:  108 on 13 and 492 DF,  p-value: <2e-16

par(mfrow = c(2, 2))
plot(fit3)

fit4 = update(fit3, ~. - age - indus)
summary(fit4)

## 
## Call:
## lm(formula = medv ~ crim + zn + chas + nox + rm + dis + rad + 
##     tax + ptratio + black + lstat, data = Boston)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.598  -2.739  -0.505   1.727  26.237 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  36.34115    5.06749    7.17  2.7e-12 ***
## crim         -0.10841    0.03278   -3.31  0.00101 ** 
## zn            0.04584    0.01352    3.39  0.00075 ***
## chas          2.71872    0.85424    3.18  0.00155 ** 
## nox         -17.37602    3.53524   -4.92  1.2e-06 ***
## rm            3.80158    0.40632    9.36  < 2e-16 ***
## dis          -1.49271    0.18573   -8.04  6.8e-15 ***
## rad           0.29961    0.06340    4.73  3.0e-06 ***
## tax          -0.01178    0.00337   -3.49  0.00052 ***
## ptratio      -0.94652    0.12907   -7.33  9.2e-13 ***
## black         0.00929    0.00267    3.47  0.00056 ***
## lstat        -0.52255    0.04742  -11.02  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.74 on 494 degrees of freedom
## Multiple R-squared:  0.741,  Adjusted R-squared:  0.735 
## F-statistic:  128 on 11 and 494 DF,  p-value: <2e-16

Nonlinear terms and Interactions

# Note: * in formula means interaction, not multiply
fit5 = lm(medv ~ lstat * age, Boston)
summary(fit5)

## 
## Call:
## lm(formula = medv ~ lstat * age, data = Boston)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -15.81  -4.04  -1.33   2.08  27.55 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 36.088536   1.469835   24.55  < 2e-16 ***
## lstat       -1.392117   0.167456   -8.31  8.8e-16 ***
## age         -0.000721   0.019879   -0.04    0.971    
## lstat:age    0.004156   0.001852    2.24    0.025 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.15 on 502 degrees of freedom
## Multiple R-squared:  0.556,  Adjusted R-squared:  0.553 
## F-statistic:  209 on 3 and 502 DF,  p-value: <2e-16

# square of lstat
fit6 = lm(medv ~ lstat + I(lstat^2), Boston)
summary(fit6)

## 
## Call:
## lm(formula = medv ~ lstat + I(lstat^2), data = Boston)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -15.28  -3.83  -0.53   2.31  25.41 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 42.86201    0.87208    49.1   <2e-16 ***
## lstat       -2.33282    0.12380   -18.8   <2e-16 ***
## I(lstat^2)   0.04355    0.00375    11.6   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.52 on 503 degrees of freedom
## Multiple R-squared:  0.641,  Adjusted R-squared:  0.639 
## F-statistic:  449 on 2 and 503 DF,  p-value: <2e-16

attach(Boston)
par(mfrow = c(1, 1))
plot(medv ~ lstat)
points(lstat, fitted(fit6), col = "red", pch = 20)
# or directly polynomical function can be used
fit7 = lm(medv ~ poly(lstat, 4))
points(lstat, fitted(fit7), col = "blue", pch = 20)

Qualitative predictors

# fix(Carseats) # you may use fix command to edit a data
names(Carseats)

##  [1] "Sales"       "CompPrice"   "Income"      "Advertising" "Population" 
##  [6] "Price"       "ShelveLoc"   "Age"         "Education"   "Urban"      
## [11] "US"

summary(Carseats)

##      Sales         CompPrice       Income       Advertising   
##  Min.   : 0.00   Min.   : 77   Min.   : 21.0   Min.   : 0.00  
##  1st Qu.: 5.39   1st Qu.:115   1st Qu.: 42.8   1st Qu.: 0.00  
##  Median : 7.49   Median :125   Median : 69.0   Median : 5.00  
##  Mean   : 7.50   Mean   :125   Mean   : 68.7   Mean   : 6.63  
##  3rd Qu.: 9.32   3rd Qu.:135   3rd Qu.: 91.0   3rd Qu.:12.00  
##  Max.   :16.27   Max.   :175   Max.   :120.0   Max.   :29.00  
##    Population      Price      ShelveLoc        Age         Education   
##  Min.   : 10   Min.   : 24   Bad   : 96   Min.   :25.0   Min.   :10.0  
##  1st Qu.:139   1st Qu.:100   Good  : 85   1st Qu.:39.8   1st Qu.:12.0  
##  Median :272   Median :117   Medium:219   Median :54.5   Median :14.0  
##  Mean   :265   Mean   :116                Mean   :53.3   Mean   :13.9  
##  3rd Qu.:398   3rd Qu.:131                3rd Qu.:66.0   3rd Qu.:16.0  
##  Max.   :509   Max.   :191                Max.   :80.0   Max.   :18.0  
##  Urban       US     
##  No :118   No :142  
##  Yes:282   Yes:258  
##                     
##                     
##                     
## 

# use all fields and interaction of Income:Advertising and Age:Price
fit1 = lm(Sales ~ . + Income:Advertising + Age:Price, Carseats)
summary(fit1)

## 
## Call:
## lm(formula = Sales ~ . + Income:Advertising + Age:Price, data = Carseats)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -2.921 -0.750  0.018  0.675  3.341 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         6.575565   1.008747    6.52  2.2e-10 ***
## CompPrice           0.092937   0.004118   22.57  < 2e-16 ***
## Income              0.010894   0.002604    4.18  3.6e-05 ***
## Advertising         0.070246   0.022609    3.11  0.00203 ** 
## Population          0.000159   0.000368    0.43  0.66533    
## Price              -0.100806   0.007440  -13.55  < 2e-16 ***
## ShelveLocGood       4.848676   0.152838   31.72  < 2e-16 ***
## ShelveLocMedium     1.953262   0.125768   15.53  < 2e-16 ***
## Age                -0.057947   0.015951   -3.63  0.00032 ***
## Education          -0.020852   0.019613   -1.06  0.28836    
## UrbanYes            0.140160   0.112402    1.25  0.21317    
## USYes              -0.157557   0.148923   -1.06  0.29073    
## Income:Advertising  0.000751   0.000278    2.70  0.00729 ** 
## Price:Age           0.000107   0.000133    0.80  0.42381    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.01 on 386 degrees of freedom
## Multiple R-squared:  0.876,  Adjusted R-squared:  0.872 
## F-statistic:  210 on 13 and 386 DF,  p-value: <2e-16

# contrasts shows how R put quantitive variables in linear regression
contrasts(Carseats$ShelveLoc)

##        Good Medium
## Bad       0      0
## Good      1      0
## Medium    0      1

---

Credit

Please note, this material is extracted from online Statistical Learning cource at Stanford University by Prof. T Hastie and Prof R. Tibshirani. It aims only for quick and future references in R and statistical learning. Please visit course page for more information and materials.

---

Linear Algebra and Modern (Markowitz) Portfolio Theory within R

In this text, i demonstrate basic linear algebra, matrix factorization (LU and QR) to solve linear equations and regression, lagrange method, basic numerical methods (bisection and Newton-Raphson) in context of modern (Markowitz, efficient frontier) portfolio theory to estimate optimum weights for given constraints within R. Firstly, lets do some basic linear equation exercises in R.

Load matrixcalc library

library("matrixcalc")

Lets assume we have a linear equation such as: Ax=b where A is square matrix and x is unknown parameters column vector and b is also column vector.

A = matrix(c(2, 2, 4, 9), 2, 2, byrow = T)
A

##      [,1] [,2]
## [1,]    2    2
## [2,]    4    9

b = c(8, 21)

To solve this equation, solve command can be used such as:

solve(A, b)

## [1] 3 1

We can also use various factorization method to resolve this linear equation.

LU Factorization

LU factorization creates elimination matrices (lower and upper) which will be used to solve linear equations.

1. The upper triangular U has the pivots on its diagonal
2. The lower triangular L has ones on its diagonal
3. L has the multipliers \( l_{ij} \) below the diagonal

lu = lu.decomposition(A)
lu

## $L
##      [,1] [,2]
## [1,]    1    0
## [2,]    2    1
## 
## $U
##      [,1] [,2]
## [1,]    2    2
## [2,]    0    5

To solve x via L (lower triangle), U (upper triangle), we have two steps:

First step is to find a column vector c which solves Lc=b equation

c = solve(lu$L, b)
c

## [1] 8 5

In the last step, Ux=c equation is solved which yields x, unknown vector.

x = solve(lu$U, c)
x

## [1] 3 1

QR Factorization

QR factorization decomposites a matrix (A) into two matrices (Q, R) , where A=QR , Q is an orthogonal matrix and R is an upper triangular matrix.

aqr = qr(A)
# orthogonal matrix
aQ = qr.Q(aqr)
aQ

##         [,1]    [,2]
## [1,] -0.4472 -0.8944
## [2,] -0.8944  0.4472

# upper triangular matrix
aR = qr.R(aqr)
aR

##        [,1]   [,2]
## [1,] -4.472 -8.944
## [2,]  0.000  2.236

Since orthogonal matrix has properties of : \( QQ^T=I \), Ax=B equation can be written as: \( Rx=Q^Tb \) where , first A decomposited into QR and then mupltiplied by \( Q^T \) .

aQT = t(aQ)
qb = aQT %*% b
x = backsolve(aR, qb)
x

##      [,1]
## [1,]    3
## [2,]    1

# all of above steps can be done with a one command too
x = qr.solve(A, b)
x

## [1] 3 1

QR factorization for Linear regression

QR factorization is also useful for linear regression fit (\( y=\beta x+\alpha \)) For example, lets try to estimate coefficients (\( \beta \) and \( \alpha \)) of linear regression between Apple and Google's montly returns.

library(quantmod)

# download price of Google and apple
getSymbols(c("GOOG", "AAPL"))

## [1] "GOOG" "AAPL"

google = c(coredata(monthlyReturn(GOOG)))
apple = c(coredata(monthlyReturn(AAPL)))
# let plot returns
plot(google, apple)

Lets apply QR factorization now to find out linear regression (apple~google)

x = cbind(1, google)
xqr = qr(x)

xQ = qr.Q(xqr, complete = T)
xR = qr.R(xqr, complete = T)
# Compute u = QTy
u = t(xQ) %*% apple
# and finally solve for $ \beta $ and $ \alpha $
backsolve(xR[1:2, 1:2], u[1:2])

## [1] 0.01636 0.66522

# lets verify the result with lm command
lm(apple ~ google)

## 
## Call:
## lm(formula = apple ~ google)
## 
## Coefficients:
## (Intercept)       google  
##      0.0164       0.6652

Linear Algebra in Portfolio Theory

Linear equations are common in modern (markowitz, efficient frontier) portfolio theory in finance to estimate optimal weights subject to some constrains, such as maximize return or minimize variance. For example:

\[ maximize: \mu^{T}w \] \[ subject to: e^{T}w=1 \] \[ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ w^T\Sigma w=\sigma^2 \]

where we have following return vector, target risk and covariance matrix :

mu = c(0.08, 0.1, 0.13, 0.15, 0.2)
mu

## [1] 0.08 0.10 0.13 0.15 0.20

s = 0.25^2  # target risk 
s

## [1] 0.0625

Sigma = matrix(c(0.0196, -0.00756, 0.01288, 0.00875, -0.0098, -0.00756, 0.0324, 
    -0.00414, -0.009, 0.00945, 0.01288, -0.00414, 0.0529, 0.020125, 0.020125, 
    0.00875, -0.009, 0.020125, 0.0625, -0.013125, -0.0098, 0.00945, 0.020125, 
    -0.013125, 0.1225), 5, 5, byrow = T)
Sigma

##          [,1]     [,2]     [,3]     [,4]     [,5]
## [1,]  0.01960 -0.00756  0.01288  0.00875 -0.00980
## [2,] -0.00756  0.03240 -0.00414 -0.00900  0.00945
## [3,]  0.01288 -0.00414  0.05290  0.02013  0.02013
## [4,]  0.00875 -0.00900  0.02013  0.06250 -0.01312
## [5,] -0.00980  0.00945  0.02013 -0.01312  0.12250

To find optimum weights, above constraint problem needs to be solved. We will use Lagrangian method for this purpose. Lagrangian is:

\[ F(w,\lambda)= \mu^{T}w + \lambda_1(e^{T}w-1)+\lambda_2(w^T\Sigma w-\sigma^2) \]

To maximize return, we can take first partial derivative of $ F'(w, \lambda)=0$ equation:

\[ G(w, \lambda) = F'(w, \lambda) = \left| \begin{array}{c} \mu+\lambda_1 e + 2\lambda_2 w \Sigma \\ e^{T}w-1 \\ w^T\Sigma w-\sigma^2 \end{array} \right| = 0 \]

To solve this quadratic equation, we can deploy numerical methods such as bisection or Newton-Raphson method. Let me explain briefly, what are these methods:

Bisection Method

If a function f(x) is continuous between [a,b] and if f(a) and f(b) have different sign, then by using following algorithm f(x)=0 can be found:

1. Increase count= count+1
2. Compute f( c ) where c =(a + b)/2 is the midpoint of [a, b]
3. If sign(f( c )) == sign(f(a)), let a=c, otherwise b=c;
4. If |b-a| > tolerance and count < totalIteration then goto 1

R implementation would be:

bisection <- function(f, a, b, tol = 0.001, iteration = 100) {
    count = 0
    while (b - a > tol & count < iteration) {
        c <- (a + b)/2
        if (sign(f(c)) == sign(f(a))) 
            a <- c else b <- c
        count = count + 1
    }
    (a + b)/2
}

This is analogous to binary search algorithm in computer science.

Newton-Raphson method

Newton-Raphson method utilises fundemantal theory of derivatives in order to find root of continuous function f(x)=0 in a recursive function, with a given starting point:

\[ x_{k+1}=x_k- \frac{f(x_k)}{f'(x_k)} \]

Compare to bisection, Newton-Raphson method converge faster, if ever converges.

Lagrange's Method & Newton-Raphson Method

Armed with Newton-Raphson method, now we can solve \( G(w, \lambda) = F'(w, \lambda) \). Algorithm to solve it:

1. Compute G(x) and G'(x)
2. Pick a starting point ($x_0$)
3. Solve the linear system : $G'(x_k)u = G(x_k)$
4. Update $x_{k+1} = x_k − u$
5. Repeat steps 3 and 4 with a number of times or $x_{k+1}- x_k$ is less than a predefined tolerance threshold.

First, we need to get partial derivative of G which will be used in in Newton-Raphson method:

\[ G'(w, \lambda) = \left| \begin{array}{ccc} 2\lambda_2 \Sigma & e & 2 \Sigma w\\ e^{T} & 0 & 0 \\ 2\Sigma w & 0 & 0 \end{array} \right| = 0 \]

Function to compute $ G(w, \lambda) $

G <- function(x, mu, Sigma, sigmaP2) {
    n <- length(mu)
    c(mu + rep(x[n + 1], n) + 2 * x[n + 2] * (Sigma %*% x[1:n]), sum(x[1:n]) - 
        1, t(x[1:5]) %*% Sigma %*% x[1:5] - sigmaP2)
}

Derivative of G

DG <- function(x, mu, Sigma, sigmaP2) {
    n <- length(mu)
    grad <- matrix(0, n + 2, n + 2)
    grad[1:n, 1:n] <- 2 * x[n + 2] * Sigma
    grad[1:n, n + 1] <- 1
    grad[1:n, n + 2] <- 2 * (Sigma %*% x[1:5])
    grad[n + 1, 1:n] <- 1
    grad[n + 2, 1:n] <- 2 * t(x[1:5]) %*% Sigma
    grad
}

Initial weights and $ \lambda $:

x = c(rep(0.5, 5), 1, 1)
x

## [1] 0.5 0.5 0.5 0.5 0.5 1.0 1.0

Lets apply Newton-Raphson iteration now:

for (i in 1:100) {
    x <- x - solve(DG(x, mu, Sigma, s), G(x, mu, Sigma, s))
}

and numerical solution:

x

## [1] -0.39550  0.09606  0.04584  0.70988  0.54372 -0.09201 -0.85715

We need to verify that $ G'(x) $ is positive definite, in order to have a maximum critical points. Lets first see what is $ G'(x) $

DG(x, mu, Sigma, sigmaP2)[1:5, 1:5]

##          [,1]      [,2]      [,3]     [,4]    [,5]
## [1,] -0.03360  0.012960 -0.022080 -0.01500  0.0168
## [2,]  0.01296 -0.055543  0.007097  0.01543 -0.0162
## [3,] -0.02208  0.007097 -0.090686 -0.03450 -0.0345
## [4,] -0.01500  0.015429 -0.034500 -0.10714  0.0225
## [5,]  0.01680 -0.016200 -0.034500  0.02250 -0.2100

Since $ G'(x) $ is symetric and square, all negative eigen value of $ G'(x) $ confirms $G'(x) is $negative definite:

eigen(DG(x, mu, Sigma, s)[1:5, 1:5])$values

## [1] -0.02024 -0.05059 -0.05806 -0.14445 -0.22364

Since all eigen values are negative, we can conclude that x[1:5] is optimum weights for the constraint and maximum return is:

t(x[1:5]) %*% mu

##        [,1]
## [1,] 0.1992

Mean reversion with Linear Regression and Bollinger Band for Spread Trading within Python

Following code demonstrates how to utilize to linear regression to estimate hedge ratio and Bollinger band for spread trading. The code can be back tested at Quantopian.com

#   Mean reversion Spread Trading  with Linear Regression
#
#   Deniz Turan, (denizstij AT gmail DOT com), 19-Jan-2014
import numpy as np
from scipy.stats import linregress

R_P = 1 # refresh period in days
W_L = 30 # window length in days
def initialize(context):
    context.y=sid(14517) # EWC
    context.x=sid(14516) # EWA
    
    
    # for long and shorting 
    context.max_notional = 1000000
    context.min_notional = -1000000.0
    # set a fixed slippage
    set_slippage(slippage.FixedSlippage(spread=0.01))
        
    context.long=False;
    context.short=False;
    
    
def handle_data(context, data):
    xpx=data[context.x].price
    ypx=data[context.y].price
    
    retVal=linearRegression(data,context)    
    # lets dont do anything if we dont have enough data yet    
    if retVal is None:
        return  None 
    
    hedgeRatio,intercept=retVal;
    spread=ypx-hedgeRatio*xpx      
    data[context.y]['spread'] = spread

    record(ypx=ypx,spread=spread,xpx=xpx)

    # find moving average 
    rVal=getMeanStd(data, context)       
    # lets dont do anything if we dont have enough data yet    
    if rVal is None:
        return   
    
    meanSpread,stdSpread = rVal
    # zScore is the number of unit
    zScore=(spread-meanSpread)/stdSpread;
    QTY=1000
    qtyX=-hedgeRatio*QTY*xpx;        
    qtyY=QTY*ypx;        

    entryZscore=1;
    exitZscore=0;

    if zScore < -entryZscore and canEnterLong(context):
        # enter long the spread
        order(context.y, qtyY)
        order(context.x, qtyX)
        context.long=True
        context.short=False    
 
    if zScore > entryZscore and canEnterShort(context):
        #  enter short the spread
        order(context.y, -qtyY)
        order(context.x, -qtyX)
        context.short=True
        context.long=False
        
    record(cash=context.portfolio.cash, stock=context.portfolio.positions_value)
    
@batch_transform(window_length=W_L, refresh_period=R_P) 
def linearRegression(datapanel, context):
    xpx = datapanel['price'][context.x]
    ypx = datapanel['price'][context.y]

    beta, intercept, r, p, stderr = linregress(ypx, xpx)
#    record(beta=beta, intercept=intercept)
    return (beta, intercept)
        
@batch_transform(window_length=W_L, refresh_period=R_P) 
def getMeanStd(datapanel, context):    

    spread = datapanel['spread'][context.y]
    meanSpread=spread.mean()
    stdSpread=spread.std()
    if meanSpread is not None and stdSpread is not None :
        return (meanSpread, stdSpread)
    else:
        return None

def canEnterLong(context):
    notional=context.portfolio.positions_value

    if notional < context.max_notional and not context.long: # and not context.short:
        return True
    else:
        return False

def canEnterShort(context):
    notional=context.portfolio.positions_value

    if notional > context.max_notional and not context.short:  #and not context.short:
        return True
    else:
        return False

Mean reversion with Kalman Filter as Dynamic Linear Regression for Spread Trading within Python

Following code demonstrates how to utilize to kalman filter to estimate hedge ratio for spread trading. The code can be back tested at Quantopian.com

#   Mean reversion with Kalman Filter as Dynamic Linear Regression
#
#   Following algorithm trades based on mean reversion logic of spread
#   between cointegrated securities  by using Kalman Filter as 
#   Dynamic Linear Regression. Kalman filter is used here to estimate hedge (beta)
#
#   Kalman Filter structure 
# 
# - measurement equation (linear regression):
#   y= beta*x+err  # err is a guassian noise 
#  
# - Prediction model:
#   beta(t) = beta(t-1) + w(t-1) # w is a guassian noise
#   Beta is here our hedge unit.
# 
# - Prediction section
#   beta_hat(t|t-1)=beta_hat(t-1|t-1)  # beta_hat is expected value of beta
#   P(t|t-1)=P(t-1|t-1) + V_w          # prediction error, which is cov(beta-beta_hat)
#   y_hat(t)=beta_hat(t|t-1)*x(t)      # measurement prediction
#   err(t)=y(t)-y_hat(t)                 # forecast error
#   Q(t)=x(t)'*P(t|t-1)*x(t) + V_e     # variance of forecast error, var(err(t))
#
# - Update section
#   K(t)=R(t|t-1)*x(t)/Q(t)                       # Kalman filter between 0 and 1
#   beta_hat(t|t)=beta_hat(t|t-1)+ K*err(t)       # State update
#   P(t|t)=P(t|t-1)(1-K*x(t))                     # State covariance update
#   
#   Deniz Turan, (denizstij AT gmail DOT com), 19-Jan-2014
#   
import numpy as np

# Initialization logic 
def initialize(context):
    context.x=sid(14517) # EWC
    context.y=sid(14516) # EWA
    
    # for long and shorting 
    context.max_notional = 1000000
    context.min_notional = -1000000.0
    # set a fixed slippage
    set_slippage(slippage.FixedSlippage(spread=0.01))
    
    # between 0 and 1 where 1 means fastes change in beta, 
    #whereas small values indicates liniar regression
    
    delta = 0.0001 
    context.Vw=delta/(1-delta)*np.eye(2);
    # default peridiction error variance
    context.Ve=0.001;

    # beta, holds slope and intersection
    context.beta=np.zeros((2,1));    
    context.postBeta=np.zeros((2,1));   # previous beta
    
    
    # covariance of error between projected beta and  beta
    # cov (beta-priorBeta) = E[(beta-priorBeta)(beta-priorBeta)']
    context.P=np.zeros((2,2));
    context.priorP=np.ones((2,2));    
    
    context.started=False;
    context.warmupPeriod=3
    context.warmupCount=0
    
    context.long=False;
    context.short=False;
     
# Will be called on every trade event for the securities specified. 
def handle_data(context, data):
    ##########################################
    # Prediction 
    ##########################################    
    if context.started:    
        # state prediction 
        context.beta=context.postBeta;
        #prior P prediction 
        context.priorP=context.P+context.Vw
    else:        
        context.started=True;
    
    
    xpx=np.mat([[1,data[context.x].price]])
    ypx=data[context.y].price
    
    # projected y
    yhat=np.dot(xpx,context.beta)[0,0]    
    # prediction error
    err=(ypx-yhat);
    # variance of err, var(err)
    Q=(np.dot(np.dot(xpx,context.priorP),xpx.T)+context.Ve)[0,0]

    # Kalman gain
    K=(np.dot(context.priorP,xpx.T)/Q)[0,0]
    
    ##########################################
    # Update section
    ##########################################    
    context.postBeta=context.beta + np.dot(K,err)

    context.warmupCount+=1
    if context.warmupPeriod > context.warmupCount:
        return
    
    #order(sid(24), 50)
    message='started: {st}, xprice: {xpx}, yprice: {ypx},\
            yhat:{yhat} beta: {b}, postBeta: {pBeta} err: {e}, Q: {Q}, K: {K}'
    message= message.format(st=context.started,xpx=xpx,ypx=ypx,\
                            yhat=yhat, b=context.beta, \
                            pBeta=context.postBeta, e=err, Q=Q, K=K)     
    log.info(message)  
   
#    record(xpx=data[context.x].price, ypx=data[context.y].price,err=err, yhat=yhat, beta=context.beta[1,0])
    ##########################################
    # Trading section
    # Spread (y-beta*x) is traded
    ##########################################    

    QTY=1000
    qtyX=-context.beta[1,0]*xpx[0,1]*QTY;        
    qtyY=ypx*QTY;        

    # similar to zscore in bollinger band 
    stdQ=np.sqrt(Q)

    if err < -stdQ and canEnterLong(context):
        # enter long the spread
        order(context.y, qtyY)
        order(context.x, qtyX)
        context.long=True
        
    if err > -stdQ and canExitLong(context):
        # exit long the spread
        order(context.y, -qtyY)
        order(context.x, -qtyX) 
        context.long=False        
 
    if err > stdQ and canEnterShort(context):
        #  enter short the spread
        order(context.y, -qtyY)
        order(context.x, -qtyX)
        context.short=True
    
    if err < stdQ and canExitShort(context):
        # exit short the spread
        order(context.y,qtyY)
        order(context.x,qtyX) 
        context.short=False
    
    record(cash=context.portfolio.cash, stock=context.portfolio.positions_value)

def canEnterLong(context):
    notional=context.portfolio.positions_value

    if notional < context.max_notional \
       and not context.long and not context.short:
        return True
    else:
        return False

def canExitLong(context):
    if context.long and not context.short:
        return True
    else:
        return False
    
def canEnterShort(context):
    notional=context.portfolio.positions_value

    if notional > context.max_notional \
       and not context.long and not context.short:
        return True
    else:
        return False

def canExitShort(context):
    if  context.short and not  context.long:
        return True
    else:
        return False

Price Spread based Mean Reversion Strategy within R and Python

Below piece of code within R and Python show how to apply basic mean reversion strategy based on price spread (also log price spread) for Gold and USD Oil ETFs.

#
# R code
#
# load price data of Gold and Usd Oil ETF 
g=read.csv("gold.csv", header=F)
o=read.csv("uso.csv", header=F)

# one month window length
wLen=22 

len=dim(g)[1]
hedgeRatio=matrix(rep(0,len),len)

# to verify if spread is stationary 
adfResP=0
# flag to enable log price
isLogPrice=0
for (t in wLen:len){
  g_w=g[(t-wLen+1):t,1]
  o_w=o[(t-wLen+1):t,1]
  
  if (isLogPrice==1){
    g_w=log(g_w)
    o_w=log(o_w)
  }
# linear regression
  reg=lm(o_w~g_w)
# get hedge ratio 
  hedgeRatio[t]=reg$coefficients[2];  
# verify if spread (residual) is stationary 
 adfRes=adf.test(reg$residuals, alternative='stationary')
# sum of p values  
  adfResP=adfResP+adfRes$p.value
}
# estimate mean p value
avgPValue=adfResP/(len-wLen)
# > 0.5261476
# as avg p value (0.5261476) indicates, actually, spread is not stationary, so strategy wont make much return. 


portf=cbind(g,o)
sportf=portf
if (isLogPrice==1){
  sportf=log(portf)
}
# estimate spread of portfolio = oil - headgeRatio*gold
spread=matrix(rowSums(cbind(-1*hedgeRatio,1)*sportf))

plot(spread[,1],type='l')

# trim N/A sections
start=wLen+1
hedgeRatio=hedgeRatio[start:len,1]
portf=portf[start:len,1:2]
spread=matrix(spread[start:len,1])

# negative Z score will be used as number of shares
# runmean and runsd are in caTools package
meanSpread=runmean(spread,wLen,endrule="constant") 
stdSpread=runsd(spread,wLen,endrule="constant")
numUnits=-(spread-meanSpread)/stdSpread #

positions=cbind(numUnits,numUnits)*cbind(-1*hedgeRatio,1)*portf

# daily profit and loss
lagPortf=lags(portf,1)[,3:4]
lagPos=lags(positions,1)[,3:4]
pnl=rowSums(lagPos*(portf-lagPortf)/lagPortf);

# return is P&L divided by gross market value of portfolio
ret=tail(pnl,-1)/rowSums(abs(lagPos))
plot(cumprod(1+ret)-1,type='l')

# annual percentage rate
APR=prod(1+ret)^(252/length(ret)) 
# > 1.032342 
sharpRatio=sqrt(252)*mean(ret)/stdev(ret)
# > 0.3713589

'''

Python code

Created on 29 Dec 2013

@author: deniz turan (denizstij@gmail.com)
'''

import numpy as np
import pandas as pd
from scipy.stats import linregress

o=pd.read_csv("uso.csv",header=0,names=["price"])
g=pd.read_csv("gold.csv",header=0,names=["price"])

len=o.price.count()
wLen=22
hedgeRatio= np.zeros((len,2))

for t in range(wLen, len):
    o_w=o.price[t-wLen:t]
    g_w=g.price[t-wLen:t]

    slope, intercept, r, p, stderr = linregress(g_w, o_w)
    hedgeRatio[t,0]=slope*-1
    hedgeRatio[t,1]=1


portf=np.vstack((g.price,o.price)).T
# spread 
spread=np.sum(np.multiply(portf,hedgeRatio),1)

# negative Z score will be used as number of shares
meanSpread=pd.rolling_mean(spread,wLen); 
stdSpread=pd.rolling_std(spread,wLen); 
numUnits=-(spread-meanSpread)/stdSpread #

#drop NaN values
start=wLen
g=g.drop(g.index[:start])
o=o.drop(o.index[:start])
hedgeRatio=hedgeRatio[start:,]
portf=portf[start:,]
spread=spread[start:,]
# number of units

numUnits=numUnits[start:,]
# position
positions=np.multiply(np.vstack((numUnits,numUnits)).T,np.multiply(portf,hedgeRatio))

# get lag 1
lagPortf=np.roll(portf,1,0);
lagPortf[0,]=lagPortf[1,];
lagPos=np.roll(positions,1,0);
lagPos[0,]=lagPos[1,];

spread=np.sum(np.multiply(portf,hedgeRatio),1)
pnl=np.sum(np.divide(np.multiply(lagPos,(portf-lagPortf)),lagPortf),1)

# return
ret=np.divide(pnl,np.sum(np.abs(lagPos),1))

APR=np.power(np.prod(1+ret),(252/float(np.size(ret,0)))) 
sharpRatio=np.sqrt(252)*float(np.mean(ret))/float(np.std(ret))
print " APR %f, sharpeRatio=%f" %( APR,sharpRatio)

Although, p-value for ADF test, APR (annual percentage rate) and sharpe ratio indicate, this strategy is not profitable, it is very basic strategy to apply.