Markov Chains in Python: a Simple Weather Model

Weather those days is so unstable and behaves strangely for the season. So, for a rainy Friday evening here is an interesting case for Python use.

Markov property is the memoryless property of a stochastic process, namely the property of future states to depend only upon the present state, not the past states. Markov chain is a process that exhibits Markov property.

Using Wikipedia example of a very simple weather model (https://en.wikipedia.org/wiki/Examples_of_Markov_chains#A_very_simple_weather_model), there is 90% probability for a sunny day tomorrow given a sunny day today and 10% probability for a rainy day tomorrow given sunny day today. Also, there is 50% probability for a sunny day tomorrow given a rainy day today and 50% probability for a rainy day tomorrow given a rainy day today. Briefly:

To solve for the weather tomorrow we need two numbers: (1) transition probabilities (as described above) and (2) initial state - what is the weather today (for instance, if today is sunny - we use (1,0), if it is rainy (0,1). Finally, we use matrix multiplication (bear in mind the difference between matrix multiplication and element-wise multiplication) to solve the problem.  The resulted state depends on previous state (changing initial state) and transition probabilities (the transition probabilities matrix remains unchanged).

Next day probability for a sunny day is 90% and probability for a rainy day is 10%. After certain number of predictions, eventually the steady state is reached and the weather prediction for the next day is independent of the initial weather. In our example at some point we reach a probability for a sunny day of 83%. and it barely changes (you can see the result below).

Here is a simple Python implementation of the weather model for a 20-day period: 

import numpy as np

transition=np.array([[0.9, 0.1],[0.5, 0.5]])

initial=np.array([1,0])

for i in range(20):

    n = 0

    counter = 0

    while counter <= n:

        counter+=1

        initial=np.dot(initial,transition)

    print (initial)

Result:

[ 0.9  0.1]

[ 0.86  0.14]

[ 0.844  0.156]

[ 0.8376  0.1624]

[ 0.83504  0.16496]

[ 0.834016  0.165984]

[ 0.8336064  0.1663936]

[ 0.83344256  0.16655744]

[ 0.83337702  0.16662298]

[ 0.83335081  0.16664919]

[ 0.83334032  0.16665968]

[ 0.83333613  0.16666387]

[ 0.83333445  0.16666555]

[ 0.83333378  0.16666622]

[ 0.83333351  0.16666649]

[ 0.8333334  0.1666666]

[ 0.83333336  0.16666664]

[ 0.83333334  0.16666666]

[ 0.83333334  0.16666666]

[ 0.83333334  0.16666666]

Portfolio Risk-Minimizing Weights Solved in Python

In a previous blogpost (Marginal and Component Value-at-Risk: A Python Example http://elenamarinova.blogspot.com/2015/05/marginal-and-component-value-at-risk.html) the Python way for calculating marginal, component and portfolio VaR was introduced. However, it left the portfolio management task a little bit incomplete as the important stage of extracting the risk-minimizing weights of the shares in the portfolio was missing.

Here is again an example of calculating VaR of a simple model portfolio with added a code for extracting weights that minimize the portfolio VaR. A change in the codes introduced in the previous post was made (we worked with annualized numbers in the previous post) to be more precise in the VaR we want, namely 99% confidence interval, 1-day VaR. (Implemented on Python 3.4)

import pandas as pd

import numpy as np

import pandas.io.data as web

import scipy.optimize as scopt

from scipy.stats import norm

Value=1e5 # $100,000

CI=0.99 # set the confidence interval

tickers =['AAPL', 'MSFT', 'YHOO']

numbers=len(tickers)

data=pd.DataFrame()

for share in tickers:

    data[share]=web.DataReader(share, data_source='yahoo', start='2011-01-01', end='2015-06-08')['Adj Close']

data.columns=tickers

ret=data/data.shift(1)-1 # calculate the simple returns

covariances=ret.cov() #gives the covariance of returns

variances=np.diag(covariances) #extracts variances of the individual shares from covariance matrix

volatility=np.sqrt(variances) #gives standard deviation

weights = np.ones(data.columns.size)/data.columns.size

weights

Out: array([ 0.33333333,  0.33333333,  0.33333333])

Pf_variance=np.dot(weights.T,np.dot(ret.cov(),weights)) #Portfolio variance

Pf_volatility=np.sqrt(Pf_variance) #Portfolio standard deviation

USDvariance=np.square(Value)*Pf_variance

USDvolatility=np.sqrt(USDvariance)

covariance_asset_portfolio=np.dot(weights.T,covariances)

covUSD=np.multiply(covariance_asset_portfolio,Value)

beta=np.divide(covariance_asset_portfolio,Pf_variance)

def VaR():

    # this code calculates Portfolio Value-at-risk (Pf_VaR) in USD-terms and Individual Value-at-risk (IndividualVaR) of shares in portfolio. 

    Pf_VaR=norm.ppf(CI)*USDvolatility

    IndividualVaR=np.multiply(volatility,Value*weights)*norm.ppf(CI)

    IndividualVaR = [ '$%.2f' % elem for elem in IndividualVaR ]

    print ('Portfolio VaR: ', '$%0.2f' %Pf_VaR)

    print ('Individual VaR: ',[[tickers[i], IndividualVaR[i]] for i in range (min(len(tickers), len(IndividualVaR)))])

VaR() # call the function

Out: Portfolio VaR:  $2830.87

Individual VaR:  [['AAPL', '$1289.91'], ['MSFT', '$1123.15'], ['YHOO', '$1493.15']]

def marginal_component_VaR():

     # this code calculates Marginal Value-at-risk in USD-terms and Component Value-at-risk of shares in portfolio.

    marginalVaR=np.divide(covUSD,USDvolatility)*norm.ppf(CI)

    componentVaR=np.multiply(weights,beta)*USDvolatility*norm.ppf(CI)

    marginalVaR = [ '%.3f' % elem for elem in marginalVaR ]

    componentVaR=[ '$%.2f' % elem for elem in componentVaR ]

    print ('Marginal VaR:', [[tickers[i], marginalVaR[i]] for i in range (min(len(tickers), len(marginalVaR)))])

    print ('Component VaR: ', [[tickers[i], componentVaR[i]] for i in range (min(len(tickers), len(componentVaR)))])

marginal_component_VaR()# call the function

Out: Marginal VaR: [['AAPL', '0.027'], ['MSFT', '0.023'], ['YHOO', '0.034']]

Component VaR:  [['AAPL', '$908.03'], ['MSFT', '$783.18'], ['YHOO', '$1139.65']]

# Risk-minimizing approach; first, we set the objective function and then we minimize it by changing weights

def VaR(weights):

    return norm.ppf(CI)*np.sqrt(np.square(Value)*(np.dot(weights.T,np.dot(ret.cov(),weights))))

def minVaR(VaR):
    nstocks = data.columns.size
    bounds = [(0.,1.) for i in np.arange(nstocks)] #non-negative weights 
    constraints = ({'type': 'eq',
                        'fun': lambda weights: np.sum(weights) - 1}) #sum of weights to equal 100%
    results = scopt.minimize(VaR, weights,
                                 method='SLSQP',
                                 constraints = constraints,
                                 bounds = bounds)
    result_weights=np.round(results.x, 4)
    result_weights= [ '%.8f' % elem for elem in result_weights ]
    new_VaR=VaR(results['x'])
             
    print ('Optimal VaR: ', new_VaR)
    print ('Optimal Weights: ', [[tickers[i], result_weights[i]] for i in range (len(tickers))])

 minVaR(VaR) # call the function

Out: Optimal VaR:  2756.04924565

Optimal Weights:  [['AAPL', '0.328'], ['MSFT', '0.475'], ['YHOO', '0.197']]

Moving Average Crossover in Python

Developing a successful algorithmic trade strategy is a challenging task, and not just because it is difficult to predict price movements based on past behavior (this is close to what we all have been taught and probably learned in our day-to-day lives). The underlying idea behind algo trading is to exploit market patterns. But market patterns do change, so fast reaction is crucial.

I was trying recently to work on some algo trade codes on Python (Python 2.7) and I ended up for a simple moving average crossover code (15-day moving average and 50-day moving average). The example is for Microsoft shares (ticker: MSFT). Moving average crossover signal is: (i) buy when 15-day moving average crossed above the 50-day moving average and (ii) sell when the 15-day moving average crossed below the 50-day moving average.

Inspired by Aaron Ash’s post “Fun With Python” (http://aaronash.com.au/posts/fun_with_python/) here is what I get:

import pandas as pd
from pandas import *
import StringIO
import urllib
import datetime
import matplotlib.pyplot as plt

%matplotlib inline

def get_data(ticker):

    base_url = 'http://ichart.finance.yahoo.com/table.csv?'

    search_query = 's={}'.format(ticker)

    search_url = '{}{}'.format(base_url, search_query)

    read_url=urllib.urlopen(search_url).read()

    return read_url

# for Python 3.4: read_url=urllib.request.urlopen(search_url).read(); for Python 2.7: read_url=urllib.urlopen(search_url).read()

if __name__ == '__main__':

    ticker = 'MSFT' #select ticker

data=get_data(ticker)

lines = data.split('\n')

lines.reverse()

lines[0] = lines[-1]

cleaned_data = "\n".join(lines[0:-1])

data = read_csv(StringIO.StringIO(cleaned_data))

data.tail() # to see how the data looks like

# calculate the 15-day and 50-day simple moving average (MA)

data['MA_15'] = rolling_mean(data.Close, 15)
data['MA_50'] = rolling_mean(data.Close, 50)
previous_15 = data['MA_15'].shift(1)
previous_50 = data['MA_50'].shift(1)

crossing = (((data['MA_15'] <= data['MA_50']) & (previous_15 >= previous_50))
| ((data['MA_15'] >= data['MA_50']) & (previous_15 <= previous_50)))
crossing_dates = data.loc[crossing, 'Date']
print(crossing_dates) #see where are the crossing points, triggering action

final=DataFrame(data, columns=['Date', 'Close', 'MA_15', 'MA_50'])

final['Date'] = pd.to_datetime(final['Date'])

for_chart=final[final['Date'] >= datetime.date(2013,1,1)]  #to see the results since Jan 1, 2013

# plot the results since Jan 1, 2013

for_chart.plot('Date', 'Close')

for_chart.plot('Date', 'MA_15')
for_chart.plot('Date', 'MA_50')
plt.legend(['Close', 'MA_15', 'MA_50'], loc='upper left')

To see the crossover dates triggers (buy or sell):

crossing_dates=pd.DataFrame(crossing_dates, columns=['Date'])

final=pd.DataFrame(data, columns=['Date', 'Close', 'MA_15', 'MA_50'])

result=crossing_dates.merge(final, on='Date')

result

result['Strategy'] = 0

condition = result['MA_15'] > result['MA_50']

result.loc[condition, 'Strategy'] = 'Buy'

result.loc[~condition, 'Strategy'] = 'Sell'

result