"""
Created on Sat Jan-03-2015
@author: Deniz Turan (http://denizstij.blogspot.co.uk/)
"""
import numpy as np
import statsmodels.api as sm
from scipy.interpolate import interp1d
def kpssTest(x, regression="LEVEL",lshort = True):
"""
KPSS Test for Stationarity
Computes the Kwiatkowski-Phillips-Schmidt-Shin (KPSS) test for the null hypothesis that x is level or trend stationary.
Parameters
----------
x : array_like, 1d
data series
regression : str {'LEVEL','TREND'}
Indicates the null hypothesis and must be one of "Level" (default) or "Trend".
lshort : bool
a logical indicating whether the short or long version of the truncation lag parameter is used.
Returns
-------
stat : float
Test statistic
pvalue : float
the p-value of the test.
usedlag : int
Number of lags used.
Notes
-----
Based on kpss.test function of tseries libraries in R.
To estimate sigma^2 the Newey-West estimator is used. If lshort is TRUE, then the truncation lag parameter is set to trunc(3*sqrt(n)/13), otherwise trunc(10*sqrt(n)/14) is used. The p-values are interpolated from Table 1 of Kwiatkowski et al. (1992). If the computed statistic is outside the table of critical values, then a warning message is generated.
Missing values are not handled.
References
----------
D. Kwiatkowski, P. C. B. Phillips, P. Schmidt, and Y. Shin (1992): Testing the Null Hypothesis of Stationarity against the Alternative of a Unit Root. Journal of Econometrics 54, 159--178.
Examples
--------
x=numpy.random.randn(1000) # is level stationary
kpssTest(x)
y=numpy.cumsum(x) # has unit root
kpssTest(y)
z=x+0.3*arange(1,len(x)+1) # is trend stationary
kpssTest(z,"TREND")
"""
x = np.asarray(x,float)
if len(x.shape)>1:
raise ValueError("x is not an array or univariate time series")
if regression not in ["LEVEL", "TREND"]:
raise ValueError(("regression option %s not understood") % regression)
n = x.shape[0]
if regression=="TREND":
t=range(1,n+1)
t=sm.add_constant(t)
res=sm.OLS(x,t).fit()
e=res.resid
table=[0.216, 0.176, 0.146, 0.119]
else:
t=np.ones(n)
res=sm.OLS(x,t).fit()
e=res.resid
table=[0.739, 0.574, 0.463, 0.347]
tablep=[0.01, 0.025, 0.05, 0.10]
s=np.cumsum(e)
eta=np.sum(np.power(s,2))/(np.power(n,2))
s2 = np.sum(np.power(e,2))/n
if lshort:
l=np.trunc(3*np.sqrt(n)/13)
else:
l=np.trunc(10*np.sqrt(n)/14)
usedlag =int(l)
s2=R_pp_sum(e,len(e),usedlag ,s2)
stat=eta/s2
pvalue , msg=approx(table, tablep, stat)
print "KPSS Test for ",regression," Stationarity\n"
print ("KPSS %s=%f" % (regression, stat))
print ("Truncation lag parameter=%d"% usedlag )
print ("p-value=%f"%pvalue )
if msg is not None:
print "\nWarning:",msg
return ( stat,pvalue , usedlag )
def R_pp_sum (u, n, l, s):
tmp1 = 0.0
for i in range(1,l+1):
tmp2 = 0.0
for j in range(i,n):
tmp2 += u[j]*u[j-i]
tmp2 = tmp2*(1.0-(float(i)/((float(l)+1.0))))
tmp1 = tmp1+tmp2
tmp1 = tmp1/float(n)
tmp1 = tmp1*2.0
return s + tmp1
def approx(x,y,v):
if (v>x[0]):
return (y[0],"p-value smaller than printed p-value")
if (v
Showing posts with label finance. Show all posts
Showing posts with label finance. Show all posts
Saturday, January 03, 2015
Stationarity Test with KPSS (Kwiatkowski–Phillips–Schmidt–Shin) in Python
In addition to augmented Dickey–Fuller test (ADF), KPSS (Kwiatkowski–Phillips–Schmidt–Shin) is widely used to verify stationarity of a signal. Python statsmodels contains ADF test, but I could not find any implementation of KPSS in python. Therefore, I converted R implementation of kpss.test in tseries library to python as follows (Click here to download the source code) :
Monday, February 24, 2014
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.
- The upper triangular U has the pivots on its diagonal
- The lower triangular L has ones on its diagonal
- 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:
- Increase count= count+1
- Compute f( c ) where c =(a + b)/2 is the midpoint of [a, b]
- If sign(f( c )) == sign(f(a)), let a=c, otherwise b=c;
- 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:
- Compute G(x) and G'(x)
- Pick a starting point ($x_0$)
- Solve the linear system : $G'(x_k)u = G(x_k)$
- Update $x_{k+1} = x_k − u$
- Repeat steps 3 and 4 with a number of times or $x_{k+1}- x_k$ is less than a predefined tolerance threshold.
\[ 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
Sunday, January 19, 2014
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
Sunday, December 29, 2013
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.
Monday, November 18, 2013
Simple Passive Momentum Trading with Bollinger Band
Below, you can see a simple trading algorithm based on momentum and bollinger band on Quantopian.com
Click here to run algorithm on Quantopian.com.
# Simple Passive Momentum Trading with Bollinger Band
import numpy as np
import statsmodels.api as stat
import statsmodels.tsa.stattools as ts
# globals for batch transform decorator
R_P = 1 # refresh period in days
W_L = 30 # window length in days
lookback=22
def initialize(context):
context.stock = sid(24) # Apple (ignoring look-ahead bias)
# for long and shorting
context.max_notional = 1000000
context.min_notional = -1000000.0
# set a fixed slippage
set_slippage(slippage.FixedSlippage(spread=0.01))
def handle_data(context, data):
# find moving average
rVal=getMeanStd(data)
# lets dont do anything if we dont have enough data yet
if rVal is None:
return
meanPrice,stdPrice = rVal
price=data[context.stock].price
notional = context.portfolio.positions[context.stock].amount * price
# Passive momentum trading where for trading signal, Z-score is estimated
h=((price-meanPrice)/stdPrice)
# Bollinger band, if price is out of 2 std of moving mean, than lets trade
if h>2 and notional < context.max_notional :
# long
order(context.stock,h*1000)
if h<-2 and notional > context.min_notional:
# short
order(context.stock,h*1000)
@batch_transform(window_length=W_L, refresh_period=R_P)
def getMeanStd(datapanel):
prices = datapanel['price']
meanPrice=prices.mean()
stdPrice=prices.std()
if meanPrice is not None and stdPrice is not None :
return (meanPrice, stdPrice)
else:
return None
Screen shot of the back testing result is:
Sunday, November 10, 2013
Cointegration Tests (ADF and Johansen) within R
In pair trading, in addition to correlation, cointegration can be very useful tool to determine which securities to be picked up. In this text, i demonstrate two approaches to test cointegration between two financial time series, two ETFs (EWA and EWC), within R. As you see below chart, in last 5 years, there is a cointegration between EWA and EWC.
- Not sure which security is dependent or independent - Can not test multiple instruments Johansen test addresses these points.
ADF Test
In this test, we use linear regression to estimate spread between two securities and then ACF to test if spread is stationary, which in a way also test of cointegration for two securities.>library("quantmod") # To get data for symbols
> library("fUnitRoots") # Unit test
## Lets get first data for EWA and EWC from yahoo finance and extract adjusted close prices
>getSymbols("EWA")
>getSymbols("EWC")
>ewaAdj=unclass(EWA$EWA.Adjusted)
>ewcAdj=unclass(EWC$EWC.Adjusted)
## Now lets do linear regression where we assume drift is zero. Since we are not sure which security is dependent and independent, we need to apply following for both case
## EWC is dependent here
> reg=lm (ewcAdj~ewaAdj+0)
## And now lets use adf test on spread (which is actually residuals of regression above step)
> adfTest(reg$residuals, type="nc")
Title:
Augmented Dickey-Fuller Test
Test Results:
PARAMETER:
Lag Order: 1
STATISTIC:
Dickey-Fuller: -1.8082
P VALUE:
0.07148
## EWA is dependent here this time
> reg=lm (ewaAdj~ewcAdj+0)
> adfTest(reg$residuals, type="nc")
Title:
Augmented Dickey-Fuller Test
Test Results:
PARAMETER:
Lag Order: 1
STATISTIC:
Dickey-Fuller: -1.7656
P VALUE:
0.07793
We use most negative Dickey-Fuller value (-1.8082 and -1.7656) to choice which regression formula to use. Based on that, We choice EWC is dependent. Within 90% confidence level (p-value is 7%), we can reject null hypothesis (unit root), so we can assume spread (residual) is stationary, therefore there is a cointegration. Below coded a function for this purpose:cointegrationTestLM_ADF <-function(A, B, startDate) {
cat("Processing stock:",A ," and ", B, " start date:",startDate)
aData=getSymbols(A,from=startDate,auto.assign = FALSE)
aAdj=unclass(aData[,6])
bData=getSymbols(B,from=startDate,auto.assign = FALSE)
bAdj=unclass(bData[,6])
lenA=length(aAdj)
lenB=length(bAdj)
N= min(lenA,lenB)
startA=0
startB=0
if (lenA!=N || lenB!=N){
startA=lenA-N+1
startB=lenB-N+1
}
cat("\nIndex start",A,":",startA," Length ",lenA )
cat("\nIndex start",B,":",startB," Length ",lenB)
aAdj=aAdj[startA:lenA,]
bAdj=bAdj[startB:lenB,]
regA=lm(aAdj~bAdj+0)
summary(regA)
regB=lm(bAdj~aAdj+0)
summary(regB)
coA <- adfTest(regA$residuals, type="nc")
coB=adfTest(regB$residuals, type="nc")
cat("\n",A," p-value",coA@test$p.value," statistics:",coA@test$statistic)
cat("\n",B," p-value",coB@test$p.value," statistics:",coB@test$statistic)
# Lets choice most negative
if (coA@test$statistic < coB@test$statistic){
cat("\nStock ",A, " is dependent on stock ",B)
cat("\np-value",coA@test$p.value," statistics:",coA@test$statistic)
p=coA@test$p.value
s=coA@test$statistic
}else {
cat("\n Stock ",B, " is dependent on stock:",A)
cat("\n p-value",coB@test$p.value," statistics:",coB@test$statistic)
p=coB@test$p.value
s=coB@test$statistic
}
return(c(s,p))
}
How to run it:
res=cointegrationTestLM_ADF("EWA","EWC",'2007-01-01')
Processing stock: EWA and EWC start date: 2007-01-01
Index start EWA : 0 Length 1731
Index start EWC : 0 Length 1731
EWA p-value 0.0501857 statistics: -1.948774
EWC p-value 0.04719164 statistics: -1.981454
Stock EWC is dependent on stock: EWA
p-value 0.04719164 statistics: -1.981454
res
-1.98145360 0.04719164
Johansen Test
As you see above ADF approach has some drawbacks such as:- Not sure which security is dependent or independent - Can not test multiple instruments Johansen test addresses these points.
> library("urca") # For cointegration
> coRes=ca.jo(data.frame(ewaAdj,ewcAdj),type="trace",K=2,ecdet="none", spec="longrun")
> summary(coRes)
######################
# Johansen-Procedure #
######################
Test type: trace statistic , with linear trend
Eigenvalues (lambda):
[1] 0.004881986 0.001200577
Values of teststatistic and critical values of test:
test 10pct 5pct 1pct
r <= 1 | 2.07 6.50 8.18 11.65
r = 0 | 10.51 15.66 17.95 23.52
Eigenvectors, normalised to first column:
(These are the cointegration relations)
EWA.Adjusted.l2 EWC.Adjusted.l2
EWA.Adjusted.l2 1.000000 1.0000000
EWC.Adjusted.l2 -1.253545 -0.3702406
Weights W:
(This is the loading matrix)
EWA.Adjusted.l2 EWC.Adjusted.l2
EWA.Adjusted.d 0.007172485 -0.003894786
EWC.Adjusted.d 0.011970316 -0.001504604
Johansen test estimates the rank (r) of given matrix of time series with confidence level. In our example we have two time series, therefore Johansen tests null hypothesis of r=0 < (no cointegration at all), r<1 (till n-1, where n=2 in our example). As in example above, if r<=1 test value (2.07) was greater than a confidence level's value (say 10%: 6.50), we would assume there is a cointegration of r time series (in this case r<=1). But as you see, none of our test values are greater than than critical values at r<0 and r<=1, therefore there is no cointegration. This is opposite of ADF result we found above. Based on my some research, i've found that Johansen test can be misleading in some extreme case (see that discussion for more info). Once a cointegration is established, eigenvector (normalized first column) would be used as weight for a portfolio.
In addition above methods, KPSS(Kwiatkowski–Phillips–Schmidt–Shin)can be also used to test stationarity.
Sunday, November 03, 2013
Stationary Tests : Augmented Dickey–Fuller (ADF), Hurst Exponent, Variance Ratio (VRTest) of Time Series within R
Several trading strategies (momentum, mean reverting, ...) are based on if data is stationary or not. In this text, i demonstrate how to test it statistically.
library("quantmod") # for downloading fx data
library("pracma") # for hurst exponent
library("vrtest") # variance ratio test
library("tseries") # for adf test
library("fUnitRoots") # for adf test
## first lets fetch USD/CAD data for last 5 years
getFX("UDSCAD")
usdCad=unclass(USDCAD) # unwrap price column
# estimate log return
n=length(usdCad)
usdcadLog=log(usdCad[1:n])
## First use Augmented Dickey–Fuller Test (adf.test) to test USD/CAD is statationary
>adfTest(usdCad, lag=1)
Title:
Augmented Dickey-Fuller Test
Test Results:
PARAMETER:
Lag Order: 1
STATISTIC:
Dickey-Fuller: 0.2556
P VALUE:
0.6978
Description:
Sun Nov 03 16:47:27 2013 by user: deniz
## As you see above, null hypothesis (unit root) can not be rejected with p-value ~70 %
## So we demonstrated it is not stationary. So if there is a trend or mean reverting. Hurst exponent (H) can be used for this purpose (Note Hursy exponent relies on that random walk diffuse in proportion to square root of time.).
#Value of H can be interpreted such as:
#H=0.5:Brownian motion (Random walk)
#H<0.5:Mean reverting
#H>0.5:Trending
> hurst(usdcadLog) Hurst exponent
[1] 0.9976377
#So, USDCAD is in trending phase.
## Another way to test stationary is to use V:
> vrtest::Auto.VR(usdcadLog)
[1] 83.37723
> vrtest::Lo.Mac(usdcadLog,c(2,4,10))
$Stats
M1 M2
k=2 22.10668 15.98633
k=4 35.03888 25.46031
k=10 56.21660 41.58861
## Another way to analyse stationary condition is via linear regression in which we will try to establish if there is a link between data and diff(data(t-1))
>deltaUsdcadLog=c(0,usdcadLog[2:n]-usdcadLog[1:n-1])
> r=lm(deltaUsdcadLog ~ usdcadLog)
> summary(r)
Call:
lm(formula = deltaUsdcadLog ~ usdcadLog)
Residuals:
Min 1Q Median 3Q Max
-0.0121267 -0.0013094 -0.0000982 0.0012327 0.0103982
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -7.133e-05 1.306e-04 -0.546 0.5853
usdcadLog 8.772e-03 5.234e-03 1.676 0.0944 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.002485 on 498 degrees of freedom
Multiple R-squared: 0.005608, Adjusted R-squared: 0.003611
F-statistic: 2.808 on 1 and 498 DF, p-value: 0.0944
> r$coefficients[2]
usdcadLog
0.008771754
## Coefficient (Beta) gives clue about if there is mean reverting. If it is negative, there is a mean reverting. As you see above, it is positive, therefore as we already concluded before, it is trending. If it was negative, we would use following to find out half life of mean revertion:
>-log(2)/r$coefficients[2]
Friday, October 25, 2013
fBasics (Basic Stats) library in R
% Load the package fBasics.
> library(fBasics)
% Load the data.
% header=T means 1st row of the data file contains variable names. The default is header=F, i.e., no names.
> da=read.table("http://www.mif.vu.lt/~rlapinskas/DUOMENYS/Tsay_fts3/d-ibm3dx7008.txt",header=T)
> ibm=da[,2] % Obtain IBM simple returns
> sibm=ibm*100 % Percentage simple returns
> basicStats(sim) % Compute the summary statistics
% Turn to log returns in percentages
> libm=log(ibm+1)*100
> t.test(libm) % Test mean being zero.
Wednesday, August 19, 2009
Pricing European Options by using Binomial Model
Binomial Model is one of the simplest pricing model for European options. Below you can find an implementation of Cox-Ross-Rubinstein (CRR) approach [1] in java for pricing of put and call European options.
With following parameters, sensitivity of the binomial pricing model as function of number of time steps can be seen in figure 1. With high number of time steps, binomial model with CRR approach converges with optimal Black-Scholes formula (e.x: 30.741 for call option)
Figure 1: Option price as a function of number of time steps.
With following parameters, sensitivity of the binomial pricing model as function of number of time steps can be seen in figure 1. With high number of time steps, binomial model with CRR approach converges with optimal Black-Scholes formula (e.x: 30.741 for call option)
Figure 1: Option price as a function of number of time steps.1: package com.denizstij.finance.pricing;
2:
3: import java.util.ArrayList;
4: import java.util.List;
5:
6: /**
7: *
8: * Estimate European Option price based on Cox, Ross and Rubinstein model
9: * @author denizstij (http://denizstij.blogspot.com/)
10: *
11: */
12: public strictfp class EuropeanOptionPricingByBinomial {
13:
14: /**
15: * Estimate European Option price based on Cox, Ross and Rubinstein model
16: *
17: * @param asset Current Asset Price
18: * @param strike Exercise Price
19: * @param volatility Annual volatility
20: * @param intRate Annual interest rate
21: * @param expiry: Time to maturity (in terms of year)
22: * @param steps : Number of steps
23: * @return Put and call price of european options based on Cox, Ross and Rubinstein model
24: */
25: public List<Double> estimatePrice(double asset,
26: double strike,
27: double volatility,
28: double intRate,
29: double expiry,
30: int steps) {
31: List<Double> results = new ArrayList<Double>();
32:
33: List<Double> stockPrices = new ArrayList<Double>();
34: List<Double> callOptionPrices = new ArrayList<Double>();
35: List<Double> putOptionPrices = new ArrayList<Double>();
36:
37: double time_step = (expiry) / steps;
38: double R = Math.exp(intRate * time_step);
39: double dF = 1 / R; // discount Factor
40:
41: double u = Math.exp(volatility * Math.sqrt(time_step)); // up boundary
42: double d = 1 / u; // down boundary (Cox, Ross and Rubinstein constraint)
43: // at leaf node, price difference factor between each node
44: double uu = u * u; // (u*d)
45: double p_up = (R - d) / (u - d); // up probability
46: double p_down = 1 - p_up; // down probability
47:
48: // initiliaze stock prices
49: for (int i = 0; i <= steps; i++) {
50: stockPrices.add(i, 0.0d);
51: }
52:
53: double sDown = asset * Math.pow(d, steps);
54: stockPrices.set(0, sDown);
55:
56: // Estimate stock prices in leaf nodes
57: for (int i = 1; i <= steps; i++) {
58: double sD = uu * stockPrices.get(i - 1);
59: stockPrices.set(i, sD);
60: }
61:
62: // estimate option's intrinsic values at leaf nodes
63: for (int i = 0; i <= steps; i++) {
64: double callOptionPrice = callPayOff(stockPrices.get(i), strike);
65: callOptionPrices.add(i, callOptionPrice);
66: double putOptionPrice = putPayOff(stockPrices.get(i), strike);
67: putOptionPrices.add(i, putOptionPrice);
68: }
69:
70: // and lets estimate option prices
71: for (int i = steps; i > 0; i--) {
72: for (int j = 0; j <= i - 1; j++) {
73: double callV = dF*(p_up * callOptionPrices.get(j + 1) +
74: p_down* callOptionPrices.get(j));
75: callOptionPrices.set(j, callV);
76: double putV = dF*(p_up * putOptionPrices.get(j + 1) +
77: p_down * putOptionPrices.get(j));
78: putOptionPrices.set(j, putV);
79: }
80: }
81:
82: // first elements holds option's price
83: results.add(callOptionPrices.get(0));
84: results.add(putOptionPrices.get(0));
85: return results;
86: }
87:
88: // Pay off method for put options
89: private double putPayOff(double stockPrice, double strike) {
90: return Math.max(strike - stockPrice, 0);
91: }
92:
93: // Pay off method for call options
94: private double callPayOff(double stockPrice, double strike) {
95: return Math.max(stockPrice - strike, 0);
96: }
97:
98: public static void main(String args[]) {
99:
100: EuropeanOptionPricingByBinomial euOptionPricing = new EuropeanOptionPricingByBinomial();
101: List<Double> results = euOptionPricing.estimatePrice(
102: 230,
103: 210,
104: 0.25,
105: 0.04545,
106: 0.5, // In terms of year
107: 10);
108: Double callOptionPrice = results.get(0);
109: Double putOptionPrice = results.get(1);
110: System.out.println("call Option Price:" + callOptionPrice);
111: System.out.println("put Option Price:" + putOptionPrice);
112: }
113: }
114:
Saturday, July 25, 2009
Fooled by Randomness, Nassim Nicholas Taleb
Think a book, whose author offends, mocks and bores its readers while he pretends he is the most clever, smart, intellectual and scientific guy in the world in each page of the book. It is not only written in a chaotic structure-- chapters, sections are not related to each other or very mixed up in terms of concepts and structure, the author also uses lots of buzzwords, redundant and wrong examples.
I am talking about "Fooled by Randomness, The hidden Role of Chance in Life and in the Market", by Nassim Nicholas Taleb. Have a look at following quote from the book:
I am talking about "Fooled by Randomness, The hidden Role of Chance in Life and in the Market", by Nassim Nicholas Taleb. Have a look at following quote from the book:
"What has more value? (a) a contract that pays you $1 million if the stock market goes down 10% on any given day in the next year; (b)a contract that pays you $1 million if the stock market goes down 10% on any given day in the next year due to a terrorist act. "
What is your answer? a or b ? Taleb claims:"I expect most people to select (b)." I am not sure IQ level of people around Taleb, but i reckon, most people would go option a. (I know that, in a normal distrubuted financial world, %10 changes in stock market is so low probability -- once every 73 to 603 trillion billion years-- , but in last 80 years, that incident happened over 2 times and five sigma deviation is over 73 times, which happens once in 7000 years if finance data is normally distrubuted. More info is here).
Taleb is a trader and scholar, works in a fixed income (bonds) financial company. He has background in science (PhD) and like many other PhD graduate, he bored in academia after some time and started to work in finance sector. In his book, he claims that probability theory is not a natural or trivial concepts for many people to comprehend. I do agree with him in this claim. But he takes his claims further and try to establish a theory and life style based on randomness (rare events) . He claims that he is expert on random events and he takes advantages of these random events in his life and business.
He exemplifies his ideas with "fictitious" characters who works as traders in finance (fixed income or equity markets). These characters are generally quite extreme and opposite of each other. For example, he kicks off the book with stock and bond market traders. The stock market trader, Steve, does not have any sound education and a risk taker and get successful so quickly. On the other hand, the bond market trader, Bob, has a degree in probability and he does not take much risk in his business (to be honest, there is not much risk in fixed income market compare to stock market). Because of calculated and not risk taking style, Bob is not rich as much as Steve.
Taleb claims that the success of the stock market trader, Steve, is based on randomness, in another words just "luck". He claims that even if you put 10.000 monkeys in stock market as trader, and monkeys trade randomly, by the end of 5 years, there will be at least a rich monkey. He also claims that after some lengthy time (10 years), all of these monkeys will disappear (The clever one would run away when he has some money, but most of them lose all of their money before they are kicked out). Taleb has a point in analogy, i think. Similar to many aspect of life, randomness also has a part in stock market. But i think, his claims that without proper analysis, research and hard work in stock market, someone would get rich randomly is just ridiculous. He dismisses that the participant of stock markets are intelligent, agile and very adaptive people. Stock market has chaotic aspects, but not totally random. It responses a deterministic way to some events (for example, if a small company merges or bought by a big company, it's share will surge). And the job of traders is to predict (or to speculate) these events in order to make profit.
I found many so-called scientific or intellectual ideas of Talebs, especially in probability theory is very mixed up. For example, he ignores the main reason of filters in statistics or engineering (filtering outliners and noise). While he claims that he takes advantages of these outliners (randomness, noise) in his real life and business, he complaints about the source of the noise for example, media and journalist (He has issue with media, TV, papers too). He contradicts himself, and gives mixed messages in different part of the book.
In brief, in his book, Taleb comes cross as an arrogant, geek person who claims 'if you disagree with me, you're an idiot and I will ignore and laugh at you.' He wrote the book for just sake of writing a book, and before clearing and organising the ideas in his mind. His writing style and personality kills the some of his nice ideas. By the end of book, i felt big disappointment. It was worst ever book i read in a long time. Luckily, the current book in my hand ("The Ascent of Money: A Financial History of the World") , is making me to forget yucky taste of this book.
Subscribe to:
Posts (Atom)