AstroML for Creating Histograms in Python

Creating histograms is not that simple as it may seem! There are of course default approaches for creating a histogram, but it is nevertheless better to know what you get instead of just looking at the resulted chart. Besides, we need not only the chart per se, but also the frequency of our data linked to respective bins. It should also be noted while it is useful to customize binning there is no guarantee it yields better result than the automatically binning.

Python offers wide range of possibilities to create histograms. But with this post I would like to introduce AstroML. AstroML (http://www.astroml.org/) is a Python module for machine learning and data mining for astronomy. 

First, let’s present some theoretical background on bin width methods. First, the two most popular rules of thumb for defining bin-width, i.e. Freedman-Diaconis and Scott and second – rules that use fitness functions, i.e. Bayesian blocks and Knuth.

The bin-width (h) and number of bins (W) under Freedman-Diaconis and Scott rules are calculated as follows: 

The other two – Bayesian blocks and Knuth’s rules – are more computationally challenging as they require minimization of a cost function. Astropy (http://astropy.readthedocs.org/en/latest/index.html) gives a brief info: Knuth’s rule chooses a constant bin size which minimizes the error of the histogram’s approximation to the data, while the Bayesian Blocks uses a more flexible method which allows varying bin widths.

Example: I use daily price data for Fondul Proprietatea (ticker: FP), a Bucharest Stock Exchange listed company. My data is stored as csv-file.

!pip install astroML

import pandas as pd

import numpy as np

from astroML.plotting import hist

data=pd.read_csv('FP.csv', delimiter=';')

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

ret= ret['FP'][~np.isnan(ret['FP'])] #remove NaN in the first row of the column (since we calculate returns we get NaN in the first day) and ‘FP’  in this code indicates the name of the column

hist(ret, bins='knuth') #Knuth’s rule

hist(ret, bins='blocks') #Bayesian blocks

hist(ret, bins='freedman') #Freedman-Diaconis

hist(ret, bins='scott') #Scott

We have our return series divided into different number of bins: ranging from 11 bins (Bayesian blocks) to 46 (Knuth’s rule).

Just for the nice math we can create our own code for bins calculation: let’s do something for the Freedman-Diaconis rule:

binwidth=2*(np.percentile(ret,75)-np.percentile(ret,25))*len(ret)**(-1.0/3.0) # bin-width is calculated in accordance with formula above.

bins=np.subtract(ret.max(),ret.min())/binwidth #this is the number of bins

plt.hist(ret.values, bins=np.linspace(min(ret), max(ret), bins))

Or alternatively (with slightly different result):

plt.hist(ret.values, bins=np.arange(ret.min(), ret.max() + binwidth, binwidth))

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

Marginal and Component Value-at-Risk: A Python Example

Value-at-risk (VaR), despite its drawbacks, is a solid basis to understand the risk characteristics of the portfolio. There are many approaches to calculate VaR (historical simulation, variance-covariance, simulation). Marginal VaR is defined as the additional risk that a new position adds to the portfolio. The practical reading of the marginal VaR, as very nicely stated in Quant at risk http://www.quantatrisk.com/2015/01/18/applied-portfolio-value-at-risk-decomposition-1-marginal-and-component-var/ is: the highest the marginal VaR, the exposure to the respective asset should be reduced to lower the portfolio VaR. The component VaR shows the reduction of the portfolio value-at-risk resulting from removal of a position. The sum of component VaR of the shares in the portfolio equals the portfolio diversified VaR, while the sum of the individual VaR presents the undiversified portfolio VaR, assuming perfect correlation between assets in the portfolio.

Also it is worth noting that based on simple mathematics, the equations can be changed, i.e. there are several ways to calculate the component VaR.

We take the following steps in Python to come up with the marginal VaR of a 3-asset portfolio:

(This works for Python 2.7 for Python 3.4 -> change in print statements in the codes. Additionally, to some of the lines some formatting can be applied - for instance returns are not presented in percentages.)

import numpy as np

import pandas as pd

import pandas.io.data as web

from scipy.stats import norm

Value=1e6 # $1,000,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-05-15')['Adj Close']

data.columns=tickers

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

ret.mean()*252 #annualize the returns

covariances=ret.cov()*252 #gives the annualized 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.random.random(numbers)

weights/=np.sum(weights) # simulating random percentage of exposure of each share that sum up to 1; if we want to plug in our own weights use: weights=np.array([xx,xx,xx])

Pf_ret=np.sum(ret.mean()*weights)*252 #Portfolio return

Pf_variance=np.dot(weights.T,np.dot(ret.cov()*252,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 to get portfolio VaR and Individual VaR in USD-terms

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 

How to Get High Frequency Stock Prices with Python (2): Netfonds Example

Getting high frequency and also real-time stock prices data (for US and Europe stocks in particular) with Python could be very simple indeed! Just use the http://hopey.netfonds.no site. Here is a simple code for generating the url-address (depending on ticker and date), reading and saving the data. Additionally, historical daily high-frequency data is available back for several days (about 15 days). 

The  url- format for the high-frequency data is (the example uses OMV shares, traded on Vienna Stock Exchange): http://hopey.netfonds.no/tradedump.php?date=20150507&paper=E-OMVV.BTSE&csv_format=txt. While for the depth the url-format is http://hopey.netfonds.no/posdump.php?date=20150507&paper=E-OMVV.BTSE&csv_format=txt. However, it seems that depth data is not available for US shares. 

Here is the code and enjoy the data:

import urllib.request

# urllib.request for Python 3.4 and urllib for Python 2.7, respectively change the def get_data(date,ticker,exchange) code below in read_url line

def get_data(date,ticker,exchange):

    base_url = 'http://hopey.netfonds.no/tradedump.php?'

    search_query = 'date={}&paper={}.{}&csv_format=txt'.format(date,ticker,exchange)

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

    read_url=urllib.request.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 = 'AAPL' # the code-format for US stocks is using directly the ticker, for instance 'AAPL' is Apple, 'MSFT' is Micrrosoft; the format for European shares is 'E-XXXX', for instance 'E-BAYND' - Bayer, 'E-TKAV' -  Telekom Austria, 'E-LLOYL' - Lloyds Banking Group

    date=20150507 # date format is YYYYMMDD

    exchange='O' # 'O' - US shares, 'BTSE' - for European shares

data=get_data(date,ticker,exchange)

with open ("trades.txt", "wb") as f:

    f.write(data)

Generating Correlated Random Numbers: Cholesky Decomposition

When it comes to generate random numbers most of the people see this exercise as a black box. In order to add some colour to this. There are several approaches to generate correlated random numbers, of which two are mainly used, namely Cholesky decomposition and spectral (eigenvalue) decomposition. These types of decomposition are used in copula estimates.

Cholesky decomposition takes the correlation (or covariance) matrix along with randomly generated numbers and correlates them. Cholesky decomposition has two, let’s call it, versions: lower and upper triangular. The Cholesky decomposition can be done in Python via Numpy and SciPy linear algebra (linalg) libraries:

(1) np.linalg.cholesky(A) # using numply linear algebra library and

(2) scipy.linalg.cholesky(A, lower=True) # using SciPy linear algebra library with lower=True indicating we want lower triangular, if we want upper triangular: lower=False. 

I prefer to use the lower triangular matrix. I also prefer working with array instead of matrix.  

Put into practice and using Python, the things look like this (we work with 3 variables):

import numpy as np

from scipy.linalg import cholesky

from scipy.stats import norm

A=np.array([[1,0.1558, 0.2227], [0.1558,1,0.5118], [0.2227,0.5118,1]]) # set the correlation matrix

B= norm.rvs(size=(3, 10000)) # 10,000 simulations to generate 3 random normally distributed variables

Cholesky = cholesky(A, lower=True) # Cholesky decomposition - lower triangular matrix

result= np.dot(Cholesky, B) # generates correlated variables

correlation=np.corrcoef(result) #check the correlation of the simulated variables

By the way, checking if the initial matrix is positive definite is also quite easy on Python:

from scipy.linalg import eigh

eigh(A) # the first array is eigenvalues and the second is eigenvectors

In our example, the correlation matrix is positive definite, meaning that all its eigenvalues are positive.  The eigenvalues of the above correlation matrix are: 0.4832, 0.8903, 1.6265.



Cholesky decomposition can also be used in the opposite case - to uncorrelate variables that are correlated.  

How to Get High Frequency Stock Prices with Python: Prague Stock Exchange Example

From time to time high frequency stock prices data is needed even for less developed CEE stock markets. I needed high frequency prices for CEZ (a Prague Stock Exchange listed company). With this task and having Python I wrote a brief code for extracting the data (the final piece of code is saving the data):

import urllib

def create_url(ticker):

    base_url = 'http://www.pse.cz/XML/ProduktKontinualJS.aspx?'

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

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

    return search_url

def clean_url(url):

    cleaned=urllib.urlopen(url).read().split("function createOnlineChart(divName)")[0]

    for character in ['k:', 'o:', 'd:new Date', '{', '}', '\r\nvar chartDataOL =', '[', ']', '(', ')', ';', ' ']:

        if character in cleaned:

            cleaned=cleaned.replace(character,'')

    return cleaned

if __name__ == '__main__':

    ticker='4169' # '4169' CEZ, '6407' Erste Group Bank, '4171' Komercni Banka, '4174' O2, '6816' NWR, '4254' Philip Morris CR

url=create_url(ticker)

cleanurl=clean_url(url)

with open ("trades.txt", "w") as f:

    f.write(cleanurl)

Copula Application in Euro Area Credit Default Swaps

Summary 

In the research we presented copula as measure of dependence that is particularly useful for identification of tail dependence between variables. We tested different copula functions for credit default swap (CDS) changes of selected euro area countries. We reached the conclusion that the Student’s t copula describes best the dependence structure of the variables. The results support the opinion that the Gaussian copula is not a suitable tool despite its widespread use. Student’s t copula has tail dependence and hence it is more useful than Gaussian copula (no tail dependence) to simulate events like joint defaults and stock market crashes. Another interesting result of the research is that in some cases we proved that certain types of Archimedean copulas provide second best fit to data (implying that the asymmetric tail dependence is also a good fit).We use the best fit copula to model market risk of CDS. 

Download file

Principal Components Analysis of the Romanian Equity Market

Principal Components Analysis (PCA) is a dimension reduction technique that aims at identifying communalities between the initial number of variables and scale down to the number to few meaningful components that are able to explain certain amount of variance. In other words, the underlying idea behind the PCA is to extract common features of the shares’ returns, i.e. instead of looking at all shares the PCA concentrates on communalities and extracts several components that are able to explain certain amount of the variability of all data. What remains unexplained by the factors is the specific risk.  

PCA extracts the factors but does not name these factors, i.e. it is only doing a partial job. But nevertheless, it is useful to identify the specific and systematic risks of the listed shares.

There are several key features that make PCA quite a sensitive technique (but as a bottom-line this is valid to all kind of techniques that deal with identifying key meaningful factors):

(1)    Identifying the number of components – a number of stopping rules apply – as the graphical analysis - scree plot (the point at which the line levels out), Kaiser’s stopping rule (only components with eigenvalue over 1.0 should be taken into consideration), etc.

(2)    Defining the options for the PCA – namely analysis of correlation or covariance and rotation method.  Generally, the covariance matrix is preferred when the variables have similar scales (for instance logarithmic stock returns) or when the variables have been transformed, while the correlation matrix is used when variables are on different scales (and correlation standardized the data). On rotation method, seems more widely accepted to use direct oblimin that allows components to be correlated (non-orthogonal solution), the alternative solution is varimax rotation (orthogonal solution – components are not correlated).

The PCA analysis of 14 Romanian listed shares for the period Sept 3, 2012 – Dec 31, 2014 (excluding several blue chips that arrived on the market after Sept 2012 – Romgaz, Electrica, Nuclearelectrica) is based on 581 daily log-return series for each company. Before conducting the PCA analysis two preliminary tests were taken into consideration – Kaiser-Meyer-Olkin Measure of Sampling Adequacy (KMO) that has a value of 0.892 (threshold value of 0.6 is considered as a minimum; closer to 1.0 the better) and Bartlett's Test of Sphericity.

The results:

(1)      Percent of variance of returns series of the 14 Romanian shares explained by the first component – 30.3% and by the second component – 7.8%; so cumulatively the first two components are able to explain as much as 38.1% of the variance of the returns. And together with the third component the % of the variance explained is 44.2%.  The first component can be interpreted as “market wide” component.

(2)    Proportion of variance of the 14 Romanian shares that can be explained by the 2 components (can be defined as the sum of the squared factor loadings) - high for SIF3, SIF1 and Transelectrica but rather low for the pharma shares – Biofarm (BIO), Antibiotice (ATB) and Fondul Proprietatea (FP).

(3)    How the companies can be combined into groups with respect to how they respond to the 2 components? Answer to this question provides the “pattern matrix”. So, in the first group we have: Transelectrica, Bucharest Stock Exchange, BRD, Banca Transilvania, Transgaz and Petrom. Almost all the large caps fall in one group (intuitive, isn’t; but this also means in case the first component is the general market we cannot rely on having sort of safe havens and really defensive stocks). Also interesting is the second group – all the SIFs (negative loadings to the second component). The SIFs have low loadings with respect to the first component (market component) – so there is opportunity for diversification between the 2 group of stocks identified. The correlation between the 2 components is -0.519.

This publication is for information purposes only and should not be construed as a solicitation or an offer to buy or sell any securities.