import math
import numpy as np
import pandas as pd
import scipy as sp
from statsmodels.tsa.vector_ar.var_model import VAR
from sklearn.covariance import OAS, EmpiricalCovariance
from statsmodels.tsa.stattools import adfuller
import statsmodels.api as smapi
def initialize(context):
set_commission(commission.PerShare(cost=0, min_trade_cost=0))
set_slippage(slippage.FixedSlippage(spread=0))
context.days = 0
context.hasPosition = False
context.stocks = None
context.positionStocks = None
context.mean = 0
context.std = 0
context.posSign = 0
context.signal = 0
context.weights = None
schedule_function(trade, date_rules.every_day(), time_rules.market_open(minutes=12))
def before_trading_start(context, data):
if context.hasPosition:
context.days += 1
return
fundamental_df = get_fundamentals(query(fundamentals.valuation.market_cap)
.filter(fundamentals.asset_classification.morningstar_sector_code == 309)
.filter(fundamentals.valuation.market_cap != None)
.filter(fundamentals.company_reference.primary_exchange_id != "OTCPK")
.filter(fundamentals.share_class_reference.security_type == 'ST00000001')
.filter(~fundamentals.share_class_reference.symbol.contains('_WI'))
.filter(fundamentals.share_class_reference.is_primary_share == True)
.filter(fundamentals.share_class_reference.is_depositary_receipt == False)
.order_by(fundamentals.valuation.market_cap.desc())).T
context.stocks = fundamental_df[0:15].index
def handle_data(context,data):
monitorPosition(context, data)
def trade(context, data):
if context.hasPosition:
return
prices = data.history(context.stocks, "price", 200,"1d")
prices = prices.dropna(axis=1)
mu = []
sig = []
BlockSize = np.shape(prices.values)[1]
insample = prices.values
for i in range(0, np.shape(insample)[0]-BlockSize, BlockSize):
x = insample[i:i+BlockSize, :]
mu.append(np.mean(x, axis=0))
sig.append(EmpiricalCovariance().fit(x).covariance_)
mu = np.asarray(mu)
sig = np.asarray(sig)
Mu = np.mean(mu, axis=0)
N = np.shape(sig)[0]
Sig = np.zeros(np.shape(sig[0, :, :]))
for i in range(0, N):
Sig += sig[i]
Sig /= N * 1.
invSig = np.linalg.inv(Sig)
S = np.zeros(np.shape(Sig))
for i in range(0, N):
S = S + np.outer(mu[i], mu[i]) + 2 * np.dot(sig[i],np.dot(invSig, sig[i]))
S = S / N - np.outer(Mu,Mu) - 2 * Sig;
phi,lambdas = sp.linalg.eig(S, Sig)
idx = phi.argsort()
stat = lambdas[:, idx[0]]
pvalue = np.dot(prices.values, stat)
if adfuller(pvalue)[1] > 0.05:
return
context.mean, context.std = fitOU(pvalue)
signal = (pvalue[-1] - context.mean) / context.std
abspval = np.dot(prices.values[-1, :], np.abs(stat))
if abs(signal) > 0.25:
context.signal = signal
print signal
openPos(-np.sign(signal), prices, stat, abspval, data, context)
def openPos(side, prices, W, abspval, data, context):
context.positionStocks = []
wsum = 0
for i, sid in enumerate(prices):
val = W[i] * data.current(sid,"price") / abspval * context.portfolio.portfolio_value * side
order_target_value(sid, val)
wsum += val
context.positionStocks.append(sid)
context.days = 0
context.posSign = side
context.weights = W
context.hasPosition = True
def monitorPosition(context, data):
if not context.hasPosition:
return
prices = []
for sid in context.positionStocks:
prices.append(data.current(sid,"price"))
pvalue = np.dot(prices, context.weights)
signal = (pvalue - context.mean) / context.std
if context.posSign > 0 and (signal > 0 or signal < -6):
closeAll(context, data, signal)
elif context.posSign < 0 and (signal < 0 or signal > 6):
closeAll(context, data, signal)
elif context.days > 35:
closeAll(context, data, signal)
def closeAll(context, data, signal):
context.positionStocks = None
context.mean = 0
context.std = 0
context.posSign = 0
context.weights = None
context.hasPosition = False
log.info(signal)
for sid in context.portfolio.positions:
order_target(sid, 0)
def fitOU(S):
n = np.shape(S)[0] - 1
Sx = np.sum(S[:-1])
Sy = np.sum(S[1:])
Sxx = np.dot(S[:-1], S[:-1])
Sxy = np.dot(S[:-1], S[1:])
Syy = np.dot(S[1:], S[1:])
a = ( n*Sxy - Sx*Sy ) / ( n*Sxx - math.pow(Sx,2))
b = ( Sy - a*Sx ) / n
sd = math.sqrt((n*Syy - math.pow(Sy,2) - a*(n*Sxy - Sx*Sy) )/n/(n-2))
delta = 1. / 252.
mu = b/(1-a)
sigma = sd * math.sqrt( -2* math.log(a)/delta/(1-math.pow(a,2)))
return mu, sigma