We think therefore we R

GARCH estimation using maximum likelihood

2012-02-09T13:57:00.000+05:30

In my previous post I presented my findings from my finance project under the guidance of Dr Susan Thomas. The results in my paper suggested that there are macroeconomic variables, particularly the INR/USD exchange rates, that help us understand the dynamics of stock returns. Although the results that I obtained were significant at 5% level, the weight of my assertion needs some more robust check before we consider the matter closed. One possible source of discrepancy we identified was that the error terms could be heteroskedastic. Meaning that one of our assumption of classical linear regression model (CLRM) estimation viz. homoskedasticity is violated. The resultant coefficient estimated, in case of heteroskedasticity of the error terms, can have underestimated standard errors, which in turn might lead to false acceptance/rejection of our null hypothesis. There are several parametric as well as non-parametric ways to remove the effect of heteroskedasticity in the error terms on the coefficient estimates. However, the method that I adopted to correct for this effect in my model was the path breaking conditional heteroskedasticity modelling propagated by Robert Engle for which he was awarded the Nobel prize in 2003.

The idea behind auto regressive conditional heteroskedasticity (ARCH) model is quite simple and straightforward. One needs to model the mean equation (this is your regression model) and along with that there has to be a specification for the modelling the squares of the errors. Suppose we have a simple AR regression:

z_t= φ*z_t-1 + e_t(Mean equation)

Now to model the conditional variance of the error terms all we need to do is a simultaneous estimation of the squares of the error terms in the following manner:

e_t = ѵ_th_t^1/2 , where ѵ_t ~NID(0,1)

h_t= α₀+ α₁e_t-1²(Variance equation)

The above specification of the mean and the variance equations is termed as a AR(1) ARCH(1) specification. This simultaneous estimation takes into account the particular form of heteroskedasticity (ARCH(1) in this case) and estimates the 'φ' coefficient accordingly. Now the above simple specification tends to pose another problem of lag selection, what lag should be considered for the variance equation. To deal with this and several other shortcomings of the simple ARCH model, Bollerslev(1986) proposed a generalized ARCH model (GARCH). The only difference being that the variance equation now becomes:

h_t= α₀+ α₁e_t-1²+ βh_t-1

Which is nothing but a GARCH(1,1) model. The beauty of this specification is that a GARCH(1,1) model can be expressed as an ARCH(∞) model. For those who are interested in learning more about ARCH and GARCH processes and the mathematics behind them here are Dr Krishnan's notes that provide an in-depth understanding on the matter. The reason why the ARCH and GARCH models rose to such prominence was because they offered us a method of not just correcting the standard errors of the estimates (as like other robust correction methods) but also provided a functional form for modelling the variance. This idea of modelling variance was heavily adopted by financial researchers. Considering the standard deviation (or variance) of asset prices as a crude estimate of volatility, financial researchers used ARCH and GARCH models to model asset price and/or other high frequency index volatility.

Now that the idea behind the GARCH modelling is clear let us plunge right into writing the conditional maximum likelihood function for a GARCH process.

Utkarsh was the generous one who provided me with the basic structure of the codes that I then customized to solve this problem at hand. Well now let me confess that a similar result or in fact a more accurate estimates of the coefficients can be obtained using unconditional maximum likelihood estimates that are offered by any of the high priced computer packages like EVIEWS, RATS, SAS etc but my motivation (apart from learning ofcourse) was to figure out a way that could give me the flexibility of modelling any functional form of heteroskedasticity. The above codes have been presented for the estimation of a GARCH(1,1) model but I could do a simple manipulation in the definition of h_tto fit any arbitrary functional form.

P.S Sincere thanks to Utkarsh for letting me share the codes. Any feedback/suggestions are welcome.

Revisiting basic macroeconomics : Illustrations with R

2012-01-06T00:45:00.000+05:30

Prologue
After 3 semesters of studying economics at IGIDR, the basics of macroeconomics still elude me. What policy shifts what curve? What determines the slope of IS-LM and AD-AS curves? What exactly was Keynes contribution to Economics? How do all these curves move in tandem? Well, I would not be surprised if you are a graduate student in economics and face similar doubts with these basic fundamentals. So in an attempt to demystify or rather simplify these macro fundamentals for econ (as well an non econ) students Isha and I have attempted this brave move of presenting the fundamentals of macroeconomics in a concise manner. We will try to be as explicit as we can but we might end up assuming some prior knowledge of economics. In such cases, we would encourage the readers active participation in improving the quality of the post wherever he/she thinks and idea is expressed with inadequate explanation.

The story
Any theory/idea, be it in pure science or economics, emanates with a hypothesis that is supported by some assumptions. As the study on the subject progresses researchers try and play with the assumptions and see if they can make it better in terms of explaining the real world as we see it. Macroeconomics is no exception, classical macroeconomics dates back to the 19th century with major proponents being Adam Smith, David Ricardo, Malthus, John Stuart Mill, J.B Say, etc. They believed in the functioning of free markets, that is a decentralized mechanism to ensure demand meets supply would result in efficient allocation of resources. Any counter-cyclical measure by a centralized government in order to smoothen the fluctuations in the business cycles would be self defeating. The centralized government, according to these classicals, was primarily to keep the budget deficit (difference between the govt. spending and tax revenue) in check. Meaning that in times of recession, when there is less revenue being generated by taxes, the government should cut spending in order to keep the deficit in check. The reason for them to propagate this idea was the if govt. has to increase spending, in times of recession, it would lead to a rise in the interest rates (since it has to borrow from the citizens) and this increase in interest rates would crowd out private investment (which is what we don't want at times of recession). This ideology was called into question during the great depression of 1929 in the US. This is when the "Say's law", which was at the heart of classical economics, failed to hold. There was mass unemployment, which was invoulantary (as opposed to the classicals view of only voluntary unemployment), insufficient demand for the produced goods, in short "Great" was an adequate adjective to describe the depressed economic scenario during the 1930's.

The Great depression of 1929 provided researchers with the much needed opportunity to revise/update the assumption that they were playing with. J.M Keynes made a break through contribution in revising the assumptions and ideologies of the classicals. He argued that market left completely to themselves might not lead to efficient allocation of resources. This was a paradigm shift in the way policy makers and researchers thought about macroeconomics. He argued that in times of recession the govt. should boost spending in order to revive the economy, contrary to the idea that the classicals propagated. Government spending and private investment were seen as complements and not substitutes. Private spending was determined not only by the interest rates in the economy but also by the expectation of future profitability (famously called the "animal spirit"). Another set of assumptions that he relaxed from the classical framework was the full flexibility of wages and prices. Classical economists believed that the market mechanism ensured that the prices of commodities adjusted instantaneously and fully to make the supply = demand. For eg. if there are 10 apples being produced and the demand turns out for only 8 apples, classical economists argues that the price of the apples would adjust in a manner (fall in this case) so that the demand rises instantaneously and exactly by 2 apples and the equilibrium is maintained. Keynes on the other hand argued that due to rigidities in the market (sticky prices, labor unions, menu costs, etc.) the price adjustment would not be instantaneous and hence the additional demand, in the short run will have to be created by the govt. by spending more. Either it can pitch in and buy the additional 2 apples or it can expedite investment in an avenue that absorbs the labor retrenched due to insufficient demand in the apple market.

Enough of literature, now let us get our hands dirty with some basic maths and visualization of some economic concept that would help us appreciate the above ideologies. We would start by illustrating the difference between the 2 idelogies using the IS-LM framework (also called the Investment Savings/Liquidity preference Money supply). The basic national income identity finds mention in the first chapters of most of the introductory Macroeconomics textbooks hence we shall start by the same identity:

Derivation

For complete and comprehensive proofs of the above equations you can refer to a text book by William Branson or another textbook by Dornbusch and Fischer.

IS curve: The points on the IS curve represent the combinations of rate of interest (i) and output (Y) for which the goods market are in equilibrium. Meaning, at these combinations of "I" and "Y", the aggregate supply of goods equals aggregate demand for goods in the economy.

LM curve: The points on the LM curve represent the combinations of rate of interest (i) and output for which the money market is in equilibrium. Meaning, at these combinations of "I" and "Y", the aggregate demand for money equals aggregate supply of money in the economy. (Note that we have the prices as exogenously given and fixed)

R codes

Simple IS-LM framework with simulated data

Equilibrium: The point of intersection of the IS and LM curve is the combination of "I" and "Y" for which both the goods and the money market are in equilibrium.

Effects of fiscal policy (or increase in government expenditure)
An increase in the government spending (fiscal expansion) results in the rightward shift of the IS curve. This happens because the autonomous component of the aggregate demand ("A" in the above derivation) increases. Fiscal policy (with exogenous price level) leads to increase in both output and interest rates.

Effect of fiscal policy

Effect of monetary policy (change in money supply)
An increase in the money supply by the central bank (monetary policy) results in the rightward shift in the LM curve. This happens because of the increase in nominal money supply (MS), and since the prices are exogenous the entire curve shifts to the right. Monetary policy (with exogenous prices) leads to fall in interest rates and rise in the output.

Effect of monetary policy

Mixture of monetary and fiscal policy
Any desired level of output and interest rates can be achieved by using both these policies in tandem, as illustrated below.

Effect of both fiscal and monetary policy

Aggregate demand curve (AD curve)
Now if we introduce prices in the picture we can trace out different combinations of "P"(prices) and "Y" for which both the goods and the money market are in equilibrium. In order to achieve this we have endogenised the computation of prices (P) (by taking ms = (nominal money)/(prices)) This represents an important relationship between aggregate price level and output for the economy. The downward sloping AD curve should not be confused with the demand curve for a good (as in microeconomics). Although both are downward sloping, the reasons for the negative relation between the prices and output demanded are different in the 2 cases. In microeconomics if the price of a good increases less of it is demanded, ceteris paribus. However in the case of AD curve, this negative relation is established by the interplay between the goods and the money market that ensures that markets clear.

This completes the story from the demand side, as to how we arrived at the aggregate demand curve for the economy. Note that we have taken a simplified version of the equations to make our point clear and to make it useful for non-econ students too. We could introduce government taxes to make the equations more realistic (and complicated) but the math looks cleaner this way and anyways the intuition remains the same even after incorporating taxes.

The next thing that we need to do is to arrive at the aggregate supply curve (labor side) for the economy, which we shall take up in the next post. Once we have presented both the ideas of aggregate demand and aggregate supply, we would be in a position to better understand the above discussion about the classical and Keynesian school of thought.

Memoization in R : Illustrative example

2012-01-04T18:32:00.000+05:30

I came across a nice problem at project euler that gave me sense of satisfaction that was unusual, I think that because I don't usually get the solutions right the first time as I did in this case. Anyhow, I shall try and decode the R codes that I used in simple English language and Mathematics.

Let me first illustrate the brute force method, that is usually the method used by novice coders like myself. The idea here is to find the largest number below 1 million that gives the maximum number of the above mentioned iterations.

So what I have done above is simply performed the iteration for each and every integer from 1 to 1 million and using a counter variable kept a track of which number gave me the largest number of iterations and recorded the corresponding number, which is what we needed in the end. The idea was straight forward the only challenge was to come up with that cute function (which we now see is not that challenging after all).

Well, now that the novice part is done lets get to what Utkarsh (my pro bro) had to say about this. My codes took ~ 701.75 seconds to run (on my Sony vaio VPCCW12EN), this was completely fine by me. Utkarsh shrugged in his usual nonchalant manner at my codes and came up with an awesome algorithm to optimize the above calculation and saving some precious time (which I think he referred to as Memoization). The idea that he worked on was that since in many cases we would already have computed the number of iterations there was no need to keep computing then again. Suppose in the example in the question we see that 13 -> 40 -> 20 -> 10 -> 5 -> 16 -> 8 -> 4 -> 2 -> 1. Now in the computation of 13 if say we already have that 10 will further iterate say 6 times we would not have to go all the way to 1. Similarly even for 10 if we know that 5 further iterates 5 times we don't need to go all the way back till 1. This would be more clear when we take a look at the codes.

The above codes, courtesy Utkarsh, took ~ 50 seconds. As it turns out I was 1,390% inefficient as compared to this optimal algorithm. I would glad to know if there is any other optimization technique (apart from using super computers) that might reduce the computational time, please share if you can find a better way of coding this in R.

Modelling returns using PCA : Evidence from Indian equity market

2011-12-26T21:15:00.004+05:30

As my finance term paper, I investigated an interesting question where I tried to identify macroeconomic variables that explain the returns on equities. Much of the debate has already taken place on this topic which has given rise to two competing theories of asset pricing viz. CAPM (capital asset pricing theory or single factor model) and APT (arbitrage pricing theory or multi-factor model). Here is a brief discussion on the two in my previous post. In this post I would like to discuss my approach to answering this question in the context of Indian stock market.

Methodology:

Companies that have been actively traded on NSE stock exchange for the past 10 years (218 companies) were selected and their daily stock returns data for these 10 years was taken from PROWESS.
Using PCA, first 10 components from the returns data of the 218 companies was extracted. More on PCA in my previous post, here.
These components were then separately regressed first on NIFTY returns (first regression)
Then these components were regressed on NIFTY returns, MIBOR rate changes, and INR/USD exchange rate changes (second regression).
The explanatory power of the 2 regressions were compared using a F-statistic. (refer to pg. 10 in the paper attached in the end of the post)

Findings and R codes:
We start with calculating the PCA of the returns on the 218 companies daily return data, then employing the 2 regressions, then comparing the 2 regressions using a F-statistic. F-stat tells us if there is any additional explanation offered when we include macroeconomic variables (viz. MIBOR, INR/USD) in our equation.

The results that I obtained pose an interesting observation. We find that the F-stat is significant at 5% for 7 out of the 10 regressions, meaning that out of the 10 regressions (each regression with a separate component) we find statistically significant addition in the explanatory power of the model after adding the macroeconomic variables. Therefore, on statistical ground I can argue that a multi factor model (APT) is preferable over a single factor model (CAPM) for modelling stock returns in the case of Indian equity market. This assertion, if holds true, can have reaching implications for asset pricing for Indian securities. Let me explain why. The principal components (that are the dependent variables in the model) are essentially the common factor across all the companies stock returns with the idiosyncratic effects discounted, so any variables that explains this common component would be the systematic risk (think why!). Now we can relate it to the debate between the CAPM and APT guys. If the CAPM guys were correct, I would obtain no additional explanation in my model after adding the macroeconomic variables i.e their assertion that the market risk (market beta) capture the entire systematic risk holds true.

The results, however, suggest that in 7 out of 10 regressions there is statistically additional explanation offered by the macroeconomic variables. Well, so we can out-rightly reject the applicability (of the much prevalent) CAPM in the case of Indian equities. Or is there something amiss? Now if I closely look at the absolute increase in the explanatory power by looking at the Adjusted-R-squared values before and after the addition of the macroeconomic variables, the absolute increase in all the cases is < 1% (refer to pg. 11 in the paper at the end of this post). Therefore, although we obtain statistical efficiency after the addition of the variables, the economic efficiency (intuition) is called to question. Is it worth while to complicate our model with additional macroeconomic variables, when we can simply have the market rate used as a reasonable proxy for all the variables? And all this just to prove a point that we have macro-variables that can provide 0.5% additional explanation in our model? This takes us back to the eternal debate of statistical vs economic efficiency, what is more important? Is the above result robust enough (on economic intuition) to question the much used, simple and powerful CAPM? Is there a threshold even in statistical efficiency to ensure economic efficiency? These are some questions that still linger on in my mind.

If we view the above result with this caveat of economic efficiency then there is reason for us to believe that a single factor model would be a preferable way to model stock returns. There are, however, evidences in the literature to suggest that multi factor (APT) is a superior way of modelling returns, but the identification of these "multi factors" remains a contentious issue among the researchers. In some desperate attempts to refute CAPM, researcher extracted principal components from a number of macroeconomic variables as the input to the PCA. This resulted in factors that had no economic intuition at all, that were then used as independent variables in explaining the returns. The APT (Arbitrage pricing theory) is a 'theory', whereas CAPM is a 'model' that approximates reality. So even if in reality there are multiple factors that give rise to the returns signals as we see them, the identification of these factors is not a trivial exercise as we have seen above. Statistically we managed to overturn the CAPM in the context of Indian equity markets but in term of economic intuition the results do not seem to be that promising. Therefore, the above exercise tells us exactly why people still stick to the evergreen CAPM as an asset pricing model.

In case you wish to replicate the exercise the data can be obtained from here: Returns_CNX_500, Nifty_returns, MIBOR, Exchange_rates.

Here is the full text of my paper. Feedback are welcome.

Movement around the mean "Stationary" OR "Unit root"

2011-12-08T02:30:00.001+05:30

The idea of modelling the time seriesof GNP, and other macroeconomic variables, data for US as a trendstationary (TS) process was brought into question by Nelson and Plosser intheir groundbreaking research paper in1982. Their research paper marked a paradigm shift in the way time-serieseconometrics was done post the 80's. The profound idea that prompted them tolook for an alternative to the prevalent TS process, was that theseries of GNP does not have any tendency to return back to a time trend following a shock. This means that following a shock (for example technologicalinnovations), the series keeps moving away from the time trend rather than returnback to it. If the series keeps moving away from the time trend, movements ofthe series would not be captured by a trend-stationary model.

This marked a radical change which transformed theidea of stationarity to include another class of processes, differencestationary (DS) processes. More on this in my previous post. But as a studentof basic time series the phenomenon of non-stationarity was not very easy forme to digest. Does it mean that if a series fluctuates around a mean, is it necessarily stationary? The answer happens to be No (now that I have completedthe course I can proudly and confidently answer that question). Accordingto the definition of stationarity, a series is stationary if any group ofconsecutive data points in the series, have the same mean. Soundsconfusing? Let me illustrate this using the example of 2 Indian macro series and some R codes.The daily 3-month MIBOR rates and the daily INR/USD exchange rates for the past10 years.

###############################

# Access the relevant files ###

###############################

mibor <- read.csv("MIBOR.csv", na.strings="#N/A")

exchange <- read.csv("Exchange_rates.csv",na.strings="#N/A")

nifty <- read.csv("Nifty_returns.csv")

#################################

## Dealing with missing values ##

#################################

## Dealing with blanks in the MIBOR rates ##

mibor[, 2] <- approx(as.Date(mibor$Dates, '%d-%b-%y'), mibor[ ,2],as.Date(mibor$Dates, '%d-%b-%y'))$y

for(k in 2:nrow(mibor)) # Calculating the%age change

{

mibor$Change1 <- diff(mibor$MIBOR) / mibor$MIBOR[-length(mibor$MIBOR)]

}

## Dealing with blanks in the exchange rates ##

exchange[, 2] <- approx(as.Date(exchange$Year,'%d-%b-%y'), exchange[,2], as.Date(exchange$Year, '%d-%b-%y'))$y

exchange$Change <- as.numeric(exchange$Change)

for(j in 2:nrow(exchange)) #Calculating the %age change

{

exchange$Change <- diff(exchange$Exchange.rates)/exchange$Exchange.rates[-length(exchange$Exchange.rates)

}

## Plotting the variables ##

png("indep_var_ns.png",width = 480, height = 480)

par(mfrow = c(2, 1))

plot(as.Date(mibor$Dates,'%d-%b-%y'),mibor$MIBOR, xlab= "Date",

ylab= "3-month MIBOR rates (%age)", type='l', col='red',

main="3-month MIBOR rates")

abline(h = 0, lty = 8,col = "gray")

plot(as.Date(exchange$Year,'%d-%b-%y'), exchange$Exchange.rates, xlab= "Date",

ylab= "IND/USD Exchange rates", type='l', col='red',

main="IND/USD Exchange rate")

abline(h = 0, lty = 8,col = "gray")

dev.off()

Eyeballing the above plots one can see that theseries do not have any trend in them, as in the series are moving more of lessabout a mean. But if we look at the MIBOR for example, the mean of the seriesis different in the period 2000-02 and different for 2003-04. This is the catchhere, which I think is quite probable to be overlooked by many. A unit root would also cause long forays away from the mean, so to conduct atest for non-stationarity we shall check if the above series has a unit root inthe auto-regressive (AR) polynomial using the ADFtest. And now that we can see that the meanis changing substantially over the time horizon, we would expect there to be aunit root in the series. Let us see what the results have to show.

> adf.test(exchange$Exchange.rates)

Dickey-Fuller= -1.9266, Lag order = 13, p-value = 0.6094 ## Cannot reject thenull of non-stationarity

alternativehypothesis: stationary

> adf.test(mibor$MIBOR)

Dickey-Fuller= -2.1925, Lag order = 13, p-value = 0.4968 ## Cannot reject thenull of non-stationarity
alternative hypothesis: stationary

> adf.test(nifty$S...P.Cnx.Nifty)

Dickey-Fuller= -11.8633, Lag order = 13, p-value = 0.00 ## Can reject the nullof non-stationarity
alternative hypothesis: stationary

So we see that the null of unit root cannot berejected for MIBOR and INR/USD, but the null is rejected for NIFTY returns. Whyits rejected for NIFTY is because the fluctuations around the mean are of avery high frequency, so even if we took 2 different time periods thestatistical difference between their means would be negligible. Thus the NIFTYreturns gives us a stationary series. MIBOR and INR/USD series are also madestationary by taking first difference of the series. The stationary plot looklike:

## Plot for the %age changes of thevariables:

png("indep_var.png", width =480, height = 480)

par(mfrow = c(3, 1))

plot(as.Date(mibor$Dates,'%d-%b-%y'),mibor$Change1, xlab= "Date",

ylab= "Change in 3-month MIBOR(%age)", type='l',col='royalblue',

main="%age change in MIBOR rates")

abline(h = 0, lty = 8, col ="gray")

plot(as.Date(nifty$Date,'%d-%b-%y'),nifty$S...P.Cnx.Nifty, xlab= "Date",

ylab= "NIFTY returns(%age)", type='l', col='royalblue',

main="NIFTY returns")

abline(h = 0, lty = 8, col ="gray")

plot(as.Date(exchange$Year, '%d-%b-%y'),exchange$Change, xlab= "Date",

ylab= "IND/USD Exchange rates change(%age)", type='l',col='royalblue',

main="IND/USD Exchange rate changes(%age)")

abline(h = 0, lty = 8, col ="gray")

dev.off()

So there are 2 takes from the exercise above (1) Series fluctuatingabout a mean need not necessarily be stationary (empirically shown) (2) 3-monthMIBOR and INR/USD exhibit unit roots in the given (10 year daily) sample for India. The first point might be a trivial statement for advancedeconometricians, but for the novice and the amateurs I think this wouldserve as a good basic exercise.

In case you wish to replicate the exercise, datacan be obtained from here: MIBOR, INR/USD, NIFTY.

Create your own Beamer template

2011-11-14T21:30:00.001+05:30

For the past couple of days, I had been searching for a tutorial that would show how to create a custom Beamer template. I found some great resources and some really great customized templates (I have listed the ones that I referred to below) but none indicated how should I go about it. There was just lines and lines of code and for someone like me, who is programmatically challenged, the task of creating my own template seemed very daunting. So I slowly started reading some of code that I found and tried to made sense of. I did eventually succeed but not entirely.

I did manage to create my own template but is relatively much simpler as compared to the ones that I saw over the internet. But this did get me started, so I am hoping that anybody who's looking for a decently documented procedure on the topic would find this post helpful. The code is relatively simple, both to understand and execute (I hope).

One of the best documents to understand Beamer is the Beamer user guide. It explains the Beamer mechanisms very well and I used it as one of the major references while preparing this template. According to the guide, Beamer presentations have five flavours of themes

Presentation themes - Every aspect of the presentation is detailed here. The colour, the font, the way the bullets look. The way the enumeration goes. The size and position of the logo etc.
Colour themes - Just of the colour details. This can be created as a separate file from the presentation theme and then later called with other presentation themes. I will explain this below. In fact, let's just enumerate the other three and jump to the explanation of how all this connects.
Font themes
Inner themes - Design the elements that are "inside" the frame like the environments, theorems, blocks etc.
Outer themes - Design the outer space of the frame like the headline, footline, sidebar etc.

To create a Beamer theme we need to specifiy four types of details-- colour, font, inner, outer. These details can be mentioned in the presentation theme itself or can be created as different files and then called in the TeX document. It is actually advisable to create them as separate files since this will allow us to use these specific themes with other presentation themes as well by just calling these in the TeX document. Additionally, it is also a more efficient way of working. For example, let's say that we like the overall feel of the Pittsburgh theme but just want to change the colours from blue to black. Then we need to just create the colour theme file and call it in the TeX document. This wouldn't be possible if we had just created one file with all the details. Let's just try to create a simple theme and see how it goes.

%%%%%%%%%%-------------------------------------------%%%%%%%%%%%%

% File created: 9 Nov, 2011

%
% This beamer style file was created just for experminetation and learning.
% The file is pretty much self-explanatory and as of 10 Nov, 2011, pretty much bug free (This is because I really don't know what "bugs" are and/or how to find them.
% Howvere, in case you find and bugs, issues, or have any suggestions/comments, please feel free to contact me at programming-r-pro-bro.blogspot.com
%
%
%
%%%%%%%%%%-------------------------------------------%%%%%%%%%%%%

% This style file is a combination of the four files required to create the Moo beamer theme.

% 1. Inner
% 2. Outer
% 3. Color
% 4. Font
% You can directly use this file instead of the using the other theme file "beamerthemeMoo-whole.sty" and calling the 4 styles in the TeX document.

%%% Defining the preamble
\mode<presentation>
\usepackage{pgfcomp-version-0-65}
\usepackage{color}

%%%%%%%%%%%%%%
%%
%% Color theme
%%
%%%%%%%%%%%%%%

%%%%%%
% We need to define a set of colours that will be assigned to various parts of the presentation. I personally believe that Beamer already does have some great themes, so if we really want to build our own, we will have go to the tiniest of details and tweak them.
%%%%%%

%%%%%%
% LaTeX, like R, has a huge set of colours to choose from and there are many ways we can access them.
% Colours in LaTeX are provided by the xcolor pacakge that loads by default. However, the xcolor package by itself does not provide enough breadth in colours.

% More colours can be called by using the "dvipsnames" options while calling beamer in the preamble of the document, e.g., \documentclass[xolor = dvipsnames]{beamer}

% Other than that we can also blend differnt colours to achive the desired colour mix or simply use the RGB codes to call that particular color.

% For blending two colours, we need to specify the two colours (Duh-uh!) and the percentage share of the two colours in the folowing format: <color1>!<percentage of color1>!<color2>
% If percentage of color1 is specified as x%, then percentage of color2 will automatically be taken as (100 - n)%.

% Why don't we try out a few options and see for ourselves.

% Here we are defining only two colours. We are primarily only going to use the first one.

\definecolor{Ftitle}{rgb}{0, 0, 0} % (rgb - 0, 0, 0) is nothing but black

% Here, we are using "rgb" in small case and this notation for decimal values of "rgb" ranging from 0 to 1

% To specify the integer values of "RGB", ranging from 0 t0 255, we need to use the "RGB" in CAPS or UPPERCASE
% Also, the first curly braces include the name that we assign to the colour combination
\definecolor{Descitem}{RGB}{0, 0, 139}

\definecolor{StdTitle}{RGB}{26, 33, 141}
\definecolor{StdBody}{RGB}{213,24,0}

\definecolor{AlTitle}{RGB}{255, 190, 190}
\definecolor{AlBody}{RGB}{213,24,0}

\definecolor{ExTitle}{RGB}{201, 217, 217}
\definecolor{ExBody}{RGB}{213,24,0}

% Another color for the background canvas using the blending option
%\definecolor{BgShade}{red!30!white}
%%% Important note: While trying this, I found out that a colour cannot be defined this way, it can only be set or called for a particular feature using the \setbeamercolor{}{} command. We will use this below to show this.

%%%%%%
% Assign colours to different constitutens of the presentation as per the requriements.
%%%%%%

% This sets the colour of the title of the presentation and titles of all the slides in the presentation to black.
\setbeamercolor{frametitle}{fg = Ftitle}
\setbeamercolor{title}{fg = Ftitle}

% In case you choose to display the Table of Contents, or the Outline slide.
\setbeamercolor{section in toc}{fg = Ftitle}
\setbeamercolor{section in toc shaded}{fg = Ftitle}

% The colour of all the items, subitems and and subsubitems are set to black.
\setbeamercolor{item}{fg = Ftitle}
\setbeamercolor{subitem}{fg = Ftitle}
\setbeamercolor{subsubitem}{fg = Ftitle}

% This sets the color for each item heading of the description environment.
\setbeamercolor{description item}{fg = Descitem}

% NOTE: Setting the color black for all the items also sets in black for other environemnts like enumerate.

% We also need to fix the colours for captions for figures and tables.
\setbeamercolor{caption}{fg = Ftitle}
\setbeamercolor{caption name}{fg = Ftitle}

% In addition, we can also change the background colour of the slides depending on pur requirement.
% \setbeamercolor{background canvas}{bg = blue!5}
% We have commented out this command because this is just for illustrative purposes and has not been used to define the background colour of the slides.

%%%%%%
% Now, there are three types of boxes in beamer:
% 1. Simple, or standard block, which can be invoked using definition or theorem
% 2. Alert block
% 3. Example block
% We will customize all these blocks based on our requirements

% Standard block

\setbeamercolor{block title}{fg = Descitem, bg = StdTitle!15!white}

\setbeamercolor{block body}{bg = StdBody!5!white}

% Alert block

\setbeamercolor{block title alerted}{bg = AlTitle}

\setbeamercolor{block body alerted}{bg = AlBody!5!white}

% Example block

\setbeamercolor{block title example}{bg = ExTitle}

\setbeamercolor{block body example}{bg = ExBody!5!white}

%%%%%%
% And one final thing, the colour of the text
\setbeamercolor{normal text}{fg = Ftitle}

%%%%%%%%%%%%%%
%%
%% Font theme
%%
%%%%%%%%%%%%%%

%%%%%%
% Here we are using default fonts
\usefonttheme{professionalfonts}

% Font for the presentation title
\setbeamerfont{title}{size = \huge}

% Font of the frame titles
\setbeamerfont{frametitle}{size = \Large}

%%%%%%%%%%%%%%
%%
%% Inner theme
%%
%%%%%%%%%%%%%%

%%%%%%
% Here we am using the rounded theme for the overall "feel" of the presentation. You can change the specific details by editing that particular option as we have done below.
\useinnertheme{rounded}

% Instead of rounded circles, we will use triangles as the indicator for items.
\setbeamertemplate{itemize items}[triangle]

% The default option for enumerate environment removes the circles around the numbers provided by the "rounded" inner theme. Just simple numbers remain.
%\setbeamertemplate{enumerate items}[default]

%%%%%%%%%%%%%%
%%
%% Outer theme
%%
%%%%%%%%%%%%%%

%%%%%%
% The outer theme takes the most amount of effort and time to customize.

%%%%%%
% In the outer theme, we will try to do the following:
% 1. Change the headline by putting a logo and a horizontal line
% 2. Change the footline and include custom information depending on our requirements
% 3. Organize the presentation title and the frame titles

%%%%%%
% Let's start with the headline. The approach that we plan to take for the headline and footline is similar.
% We will first define a new command and then include the command in the \setbeamertemplate{} option.
% I tried a couple of approaches that would make the method more simple but could not come up with one. % In case you do find a more aesthetic approach, please do send it accross. The contact information is on the top of the page.
% First the horizontal line on the top portion of the slides
% Add a horizontal line that runs from left of the slide to the right, just below the logo.
\newcommand{\LogoLine}{%
\raisebox{-12mm}[0pt][0pt]{%
\begin{pgfpicture}{0mm}{0mm}{0mm}{0mm}
\pgfsetlinewidth{0.28mm}
\color{gray}
\pgfline{\pgfpoint{-3mm}{1mm}}{\pgfpoint{10.8cm}{1mm}}
\end{pgfpicture}}}

% Include the line that we just created in the headline
\setbeamertemplate{headline}[text line]{\LogoLine}

% Now the logo. As it turns out. I could not include the logo in the headline. If I tried to, the headline kept shifting downwards.
% Acting a little smart and lazy, I just included the logo in the right sidebar and shifted it up.
\setbeamertemplate{sidebar canvas right}{
\vspace*{3pt}\hspace*{-25pt}%
{\includegraphics[height=28pt]{moo.png}}}

%%%%%
% Now that we have changed the headline, we will need to orient the frame titles in a way that the come at the right stop just above the horizontal line.
\setbeamertemplate{frametitle}{
\vspace*{4mm}\hspace*{-2mm}\insertframetitle}

%%%%%
% Like mentioned above, we will take a similar approach to customize the footline as well and inclued FAA in it.

\newcommand{\Ffootline}{%

\insertsection % The left end of the footline

\hfill

\textit{Moo} % The center

\hfill

\insertframenumber/\inserttotalframenumber} % And the right end

\setbeamertemplate{footline}{%

\usebeamerfont{structure}

\begin{beamercolorbox}[wd=\paperwidth,ht=2.25ex,dp=1ex]{title in head/foot}%

\Tiny\hspace*{4mm} \Ffootline \hspace{4mm}

\end{beamercolorbox}}

%%%%%%

% We will also remove the navigation symbols, which I personally don't find very useful
\setbeamertemplate{navigation symbols}{}
%%%%%%
% Now the toughest part--at least for me--customizing the title page.
% Putting a logo on the title page and text beside it was quite a difficult task and to be fairly honest, my code is not effecient at all. At the very least, it is plain clumsy.
% Though it does solve the purpose (I hate to use this phrase) it is neither "neat" nor "cool".
% Anyway, let's see how we went about it.
% For the title page, we needed a logo on the left, a vertical separater line, and finally a place for the title, author, date etc.

% First, let's create the line
\newcommand{\TitleLine}{%
\raisebox{-12mm}[0pt][0pt]{%
\begin{pgfpicture}{0mm}{0mm}{0mm}{0mm}
\pgfsetlinewidth{0.10mm}
\color{gray}
\pgfline{\pgfpoint{55mm}{0mm}}{\pgfpoint{55mm}{50mm}}
\end{pgfpicture}}}

% Now let's create commands for the title etc., that we can call later

% Title
\newcommand{\MyTitle}{%
\hspace*{60mm}\vspace{-25mm}
\centering \inserttitle}

% Subtitle
\newcommand{\MySubTitle}{%
\hspace*{60mm}\vspace{-25mm}
\centering \footnotesize \textit{\insertsubtitle}}

% Author
\newcommand{\MyAuthor}{
\hspace*{60mm}\vspace{-25mm}
\centering \insertauthor}

% Institute
\newcommand{\MyInstitute}{
\hspace*{60mm}\vspace{-25mm}
\centering \footnotesize \textit{\insertinstitute}}

% Date
\newcommand{\MyDate}{
\hspace*{60mm}\vspace{-25mm}
\centering \insertdate}

% We declare the image that will be used as the logo
\pgfdeclareimage[width = 0.20\paperwidth]{big}{moo.png}

% This is quite a complicated command. We basically create a "beamercolorbox" for each field and invoke the commands that we had created earlier.

\setbeamertemplate{title page}{\TitleLine

\hspace*{11mm}\vspace*{-60mm}

\begin{beamercolorbox}[wd=0.5\paperwidth,ht=0.13\paperwidth]{Title}

\pgfuseimage{big}

\end{beamercolorbox}

\begin{beamercolorbox}[wd=\paperwidth,ht=0.06\paperwidth]{Title}

\usebeamerfont{Title}%

\MyTitle

\end{beamercolorbox}

\begin{beamercolorbox}[wd=\paperwidth,ht=0.03\paperwidth]{Title}

\usebeamerfont{Title}%

\MySubTitle

\end{beamercolorbox}

\begin{beamercolorbox}[wd=\paperwidth,ht=0.06\paperwidth]{Title}

\usebeamerfont{Title}%

\MyAuthor

\end{beamercolorbox}

\begin{beamercolorbox}[wd=\paperwidth,ht=0.03\paperwidth]{Title}

\usebeamerfont{head/foot}%

\MyInstitute

\end{beamercolorbox}

\begin{beamercolorbox}[wd=\paperwidth,ht=0.07\paperwidth]{Title}

\usebeamerfont{Title}%

\MyDate

\end{beamercolorbox}}

\mode

<all>

Well, this pretty much completes the creation of Moo. I will upload the Sty file, TeX file and sample pdf for direct use and references.

Also, in case you are trying to create a more complicated theme you should definitely check out some of the resources mentioned below.

Most important - The Beamer user guide
A great repository of themes here
A great introduction to Beamer here
Some information that helped
Another custom theme here
A very good Beamer example

Hope this information helps.

The uploaded file

Style file

TeX File

PDF output

Unit root versus breaking trend: Perron's criticism

2011-11-05T02:38:00.001+05:30

I came across an ingenious simulation by Perronduring my Time-series lecture which I thought was worth sharing. The idea wasto put your model to a further test of breaking trend before accepting the nullof unit root. Let me try and illustrate this in simple language.

A non-stationary time series is one that has itsmean changing with time. In other words, if you randomly choose a bunch of values fromthe series from the middle, you would end up with different values of mean fordifferent bunches. In short there is a trend in the data which needs to beremoved to make it stationary and proceed with our analysis (its far easier towork with stationary timeseries). In order to deal with non-stationarytime-series one has to be careful about the kind of non-stationarity thatis exhibited by the variable. Two corrections for non-stationarity include fitting

(1) Trend stationary (TS) models, whichare suitable for models that have a deterministic trend and fluctuations aboutthat deterministic trend. This can be fit by a simple z_t = a + bt + e_twhere e_t~ARMA(p,q)

(2) Difference stationary (DS) models, whichare suitable for models having a stochastic trend. The DS models areappropriate for models that have a unit root in the AR polynomial. Unit root in theAR polynomial means that the trend part in the series cannot be represented bya simple linear trend with time (a + bt). And the correct representation is (1 – B)z_t= a + e_t,where e_tis i.i.d.

The asymptotic properties of the estimates, forecasts andforecast errors vary substantially between the TS and DS models. (For the onesinterested in the algebra behind this, lecture notes of Dr. Krishnan are here) Therefore it is important for us be sure thatthe model belongs to the appropriate class before we fit a TS or DS model. Thisis the reason why the clash between the two school of thoughts has bred enormousliterature and discussions on the methodology to check for unit roots. Onecould try and endlessly argue about these discussions but I want toillustrate the genius of Perron who criticized the idea of fitting a DS modelto series that could have a structural breaks. He said that you ought to takeinto account the structural break before you check for the unit roots, if youdon't do so, you might end up accepting the null of unit root, even when thetrue data generating process (DGP) is a trend stationary process. Heillustrated this using a simple, but very elegant, simulation exercise. Madhavand I, along with fine-tuning on the codes provided by Utkarsh, replicated thisexercise with R.

The steps involved are as follows:

(1) Simulate 1000 series with the DGP as:

z_t= u₁ + (u₂– u₁)DU_t+bt + e_t

where e_tare i.i.d innovations and t = 1,2,3,...100. Forsimplicity I have assumed b = 1 and u₁= 0.

(2) Assume that there is a crash at time T_b= 50 and the entire series comes down by amount u₂.

## Simulating a trend stationary model with a crash in the intercept ##

t <- c(1:100) # specifying the time

dummy <- as.numeric(ifelse(t <= 50, 0, 1)) # specifying the dummy for trend break at T = 50

z <- ts(t - 15*dummy + rnorm(100, mean = 0, sd = 3))# This is the trend stationary model with break in trend

x <- ts(t - 15*dummy) # This is just the trend line that we see in "red" in the plot below

plot(z, main = "Simulated series with break at T = 50")

lines(x, col = "red") ## Plotting a sample of the model that we have simulated

(3) For these simulations compute theautoregressive coefficient, “rho” in the regression:

z_t = u + bt + ‘rho’z_t-1+ e_t

(4) Plot the cumulativedistribution function (c.d.f) of “rho” for different values of u₂(crash).

## Now we will simulate the sample data above 1000 times and check for unit roots for each of these samples ##

# For simplicity we define a function to generate the "rho's" for each of the simulated series

sim <- function(crash) ## Function name "sim"

{

d <- ts(t - crash*dummy + rnorm(100, mean = 0, sd = 3))

## saving the simulated series in "d"

trend <- lm(d ~ t) ## remove the trend from the

simulated series

# crash in the above function refers to the value of u₂in equation 1

res <- ar(ts(trend$residuals), order=1, aic= FALSE) ##

Fit an AR(1) model to the residue obtained after

detrending the series

if(length(res$ar) < 1) 0 else res$ar ## Return the ar

coefficient of the fitted AR(1) model above.

}

## Generate "rho's" for different magnitude of crash by

simply using the sim() function defined above

temp1 <- replicate(n, sim(10))

temp2 <- replicate(n, sim(15))

temp3 <- replicate(n, sim(20))

temp4 <- replicate(n, sim(35))

## Sort the values of "rho", we do this to plot the CDF

as we will see shortly

temp1.1 <- sort(temp1)

temp2.1 <- sort(temp2)

temp3.1 <- sort(temp3)

temp4.1 <- sort(temp4)

y <- seq(from=0, to=1, length.out=n)## This is how I

define the y-axis of my CDF which are basically the

probabilities.

## Plotting all the CDF of rho for different magnitude in one plot.

plot(c(min(temp1.1), max(temp4.1)), c(0, 1), type='n', xlab = "Rho", ylab= "Probability", main = "CDF of 'Rho' for differnt magniturde of crashes")

lines(temp1.1, y, type = 'l', col = 'red')

lines(temp2.1, y, type = 'l', col = 'green')

lines(temp3.1, y, type = 'l', col = 'blue')

lines(temp4.1, y, type = 'l', col = 'black')

b <- c("10 unit crash", "15 unit crash", "20 unit crash", "35 unit crash")

legend("topleft", b , cex=0.5, col=c("red", "green", "blue", "black"), lwd=2, bty="n")

An interesting observation that we make (or ratherPerron made) is that the c.d.f of our autoregressive coefficient “rho” tendsmore towards unity with increase in the magnitude in crash. What this means isthat as the magnitude of crash increases the possibility of your accepting the (false)null of unit root increases. Why I say the false null is because I know thetrue DGP is a trend stationary one.

This idea of Perron was criticised on the groundthat he was specifying the break point (T_b) exogenously, that isfrom outside the DGP. Frankly speaking I do not understand why was this takenas a criticism. I think fixing the break point exogenously was a good way offixing it with an economic intuition and not making is a purely statisticalexercise. Some researchers (I don’t understand why) termed this (simulation) illustrationas a “data mining” exercise, and improved it by selecting the break point (T_b)endogenously (by Zivot and Andrews as mentioned in the lecture notes).

I would hate to impose my opinion here but I feelthis was a very elegant and logical way of driving home the point that the nullof unit root should be accepted for your sample if and only if your model standsthe test of extreme rigour and not otherwise, and the rigour could be imposedexogenously with economic intuition too.

P.S. Perron did a similar simulation for breakingtrend model, i.e where the slope of the model had a structural break. The codeswould be quite similar to the ones given above, in fact it would be a goodpractice if you could do the similar simulation for a breaking trend. In case you do want to try but face any issues please feel free to post/email your queries.

Criticism and discussions welcome.

Modern Portfolio Optimization Theory: The idea

2011-11-03T23:03:00.000+05:30

We were recently given a lecture (by Dr. Susan Thomas) on Harry Markowitz portfolio optimization theory, and I was really fascinating with the noble laureate's story of how he found it difficult to convince his guide about the importance of his thesis work. Little did anyone know that his thesis would get him the most respected award in academia 35 years down the lane. Let me try and illustrate what was the basic idea behind the modern portfolio theory in simple English language. The principal task of any financial advisory firm/individual is to try and solve the capital allocation problem of an agent, meaning how much money should be invested in what asset/security to get the desired rate of return with minimum risk. Lets try and visualize this problem of capital allocation from the individuals point of view. Suppose I want a 15% p.a rate of return on my investment and there is a security in the market that gives you an expected 15% annual rate of return with some risk(sigma or standard deviation) associated with that security/asset. One way for me to achieve the desired expected return is to invest in this security take the risk(sigma) associated with this security and hope that I am awarded for the risk I took and get a 15% return in the end. However, the beauty of financial markets is that since there are so many assets/securities being traded in the market, it is possible for me to diversify my investment, in simple words, park fraction of my money with different assets in a way that gives me approximately ~ 15% return but with substantially less risk. This is referred to as "diversification" in Finance parlance. The mathematics behind the results can be easily shown, one can find the maths in the lecture slides here. Let me try and illustrate this result through empirical simulation using R.

Suppose you are in a universe of only risky assets, i.e all assets that have a positive sigma. Then what Harry Markowitz did, with mathematics, was to compute the set of Expected returns that you could achieve with minimum risk (sigma). To illustrate this result I took stock returns for these 7 random companies for about 251 days, computed their respective expected returns, and the variance-co-variance matrix of the returns. With this I have all the necessary fodder to compute the E(r) and sigma for all the possible combination of weights that I can assign between these securities, meaning all possible different combinations of parking my money among these securities.

# Map your working directory using setwd()
# Read the relevant file
r1 <- read.csv("Markowitz.csv")

mu = colMeans(r1[,-7]) ## Calculate the column means and storing in "mu", I have removed the 7th column from the data frame as I have a risk free asset's return there, I will come back to this later in the post.

bigsig = cov(r1[,-7]) ## Variance-co-variance matrix

m = nrow(bigsig)-1
w = diff(c(0,sort(runif(m)), 1)); ## Assigning random weights between 0-1 to "w" which will have the dimension = no. of securities between which you have to divide your investment.

rb = sum(w*mu); ## Creating matrix "rb" which stores the E(r)

sb = sum(w*bigsig*w); ## Creating matrix "sb" which stores sigma's

N = 2000 ## Number of different combinations of "w" you want to look at

## Simulating the different combinations of weights "w's" and storing the E(r) and sigma^2

for (j in 2:N) {

w = diff(c(0,sort(runif(m)), 1));

r = sum(w*mu); rb = rbind(rb,r); 

s = sum(w*bigsig*w); sb = rbind(sb,s); ## Note this is sigma^2 (variance)

}

d = data.frame(rb, sb); ## Merge all the E(r) and sigmas in one data.frame.

d$sb = sqrt(d$sb); ## Square root the variance to get the sigmas

plot(d$sb, d$rb, ylab="E(r)", xlab="Sigma", col="blue",  xlim = c(0.5,10), ylim = c(-0.1,0.5), main = "E(r)- sigma (With risk free asset), N = 200000")

With only Risky assets

I did this simulation for 200000 different values of weights too, for those who are curious the plot can be obtained from here.

If we close the envelope of the above mentioned points with a smooth curve (from the left hand side) we get the efficient portfolio frontier (EFF). So what the graph suggests is that one would want to choose the weights (allocations) that lie, more specifically, on the above half of the curve. It makes intuitive sense too, why would one take on more risk to get a lesser return (which is implied by the lower half of the curve). Also another observation that we can make from the above plot is that for every additional unit of returns we will have to take on an increasingly higher risk (illustrated by the upper half of the plot). This risk-return trade off graph changes its characteristics when a risk free assets comes into the picture. Lets see what happens is we incorporate the risk free asset in our data frame. Note that a risk free assets could be thought of as a government security, whose risk (sigma) is 0, hence even its co-variances with the other risky assets would be 0. The only addition to the above codes would be to replace r1[,-7] with just r1 that is include the risk free asset too.

With risk-free asset

Simulation with N = 200000 here.

Look what happened to our risk-return trade off relationship, its linear!! So its a straight line instead of a smooth curve. Meaning for every additional unit of return you will have to take a fixed additional amount of risk, which will be given by the slope of the above half of the graph. This result, which sounds to be a trivial finding, can be extended to the entire universe of assets and one can compute the optimal portfolio that will give you a desired rate of return with minimum sigma. However the implementation of this methodology in the actual world is the real challenge. Computation of E(r), sigma, getting the up to date data for all the securities or assets in the universe are some of the challenges that one would face. But nevertheless it is inspiring to see how Harry Markowitz stuck to his idea and triumphed in the end with this path breaking research that shaped the financial advisory industry as we see it now.

P.S In case you want to replicate the exercise the data can be downloaded from here.

Modelling with R: part 5

2011-10-30T19:27:00.001+05:30

In our exercise of learning modelling in R, we have till now succeeded in doing the following:

Specifically, we created a decision tree using the rpart package. The decision tree was built with no priors and gave an AUC value of 0.73475 as opposed to 0.80493 given by the logistic regression. So here the logistic regression outperforms the recursive partitioning methodology of the rpart package. In this post we will try to see how we can improve over the current performance level of our decision tree. As a way to gauge our performance in the future, we can consider the AUC of the logistic regression as a minimum benchmark over which we need to improve. So let's get started.

## 2.2.2: Decision tree with priors ##

t1.prior <- rpart(default ~

amt.fac +

age.fac +

duration +

chk_acct +

history +

purpose +

sav_acct +

employment +

install_rate +

pstatus +

other_debtor +

time_resid +

property +

other_install +

housing +

other_credits +

job +

num_depend +

telephone +

foreign

, data = dev, parms = list(prior = c(0.9, 0.1)))

# Not the difference in the commands for a tree with priors and a tree without one. Here we need to specify the priors along with the formula in the rpart() function command.

Before we go further it is important to understand what "priors" exactly are. I initially did not understand the concept of priors until I read this document. You need not go through the entire document. I have made a little summary for myself that might help.

Simply put, priors are nothing but a weighing scheme. They specify the weight that we put on the overall level of the dependent variable. Technically, all trees are built with priors. In cases where priors are not specified, they are automatically taken to be proportional response rate, i.e, for each class they are proportional to the to the number of records in that class.

Take for example our data set. The dependent variable has two categories - 0 and 1. There are 300 1s and 700 0s, so the response rate is about 30% (300/(700 +300)). Then most decision trees will automatically apply priors of 30% and 70% to the respective classes. Now a simple decision tree tries to maximize the number of cases it classifies correctly and hence a lot of effort will concentrated on classifying the non-defaulters instead of the defaulters. This is because even if we incorrectly classify all the bad loans we will be wrong 30% of the time as opposed to 70% of time if we incorrectly classify the non-defaulters. As mentioned in the document, it is important to note that such a model will be, "literally accurate but practically worthless."

How can we then deal with such situations. There are three possible options here

Increase the proportion of defaulters by copying the observations so that there are 700 instances of each.

Duplicating observations, however, is not a very neat idea.
Additionally, it can be a problem if the data set is large. Adding observations will only make the other process computationally intensive and/or plain inefficient.

Randomly pick 300 non-defaulters and create a data set which has 300 instances of each.

This involves throwing out a large chuck data that contains valuable information and hence the results may not be accurate

Use a smarter decision tree algorithm that let's you specify priors.

In case you haven't already guessed we have taken the third option. If we can vary the priors based on our understanding, we can analyze data without any special handling, regardless of how skewed the distribution of the dependent variable may be.

Coming back to our code

plot(t1.prior)

# Plots the trees

text(t1.prior)

# Adds the labels to the trees.

We don't need to prune this model and can score it right away

val$t1.p.yhat <- predict(t1.prior, val, type = "prob")

We can plot the ROC curve for the tree with priors.

t1.p.scores <- prediction(val$t1.p.yhat[,2], val$default)

t1.p.perf <- performance(t1.p.scores, "tpr", "fpr")

# Plot the ROC curve

plot(t1.p.perf, col = "blue", lwd = 1)

# Add the diagonal line and the ROC curve of the logistic model, ROC curve of the tree without priors

plot(m1.perf, col = "red", lwd = 1, add = TRUE)

plot(t1.perf, col = "green", lwd = 1.5, add = TRUE)

abline(0, 1, lty = 8, col = "grey")

legend("bottomright", legend = c("tree w/o prior", "tree with prior", "logit"), col = c("green", "blue", "red"), lwd = c(1.5, 1, 1))

KS statistic

ks1.p.tree <- max(attr(t1.p.perf, "y.values")[[1]] -(attr(t1.p.perf, "x.values")[[1]]))

ks1.p.tree

AUC

t1.p.auc <- performance(t1.p.scores, "auc")

t1.p.auc

Predictability of stock returns : Using acf()

2011-10-28T01:33:00.000+05:30

In my previous post, I employed a rather crude and non-parametric approach to see if I could predict the direction of stock returns using the function runs.test(). Lets go a step further and try modelling this with a parametric econometric approach. The company that I choose for the study is INFOSYS (NSE code INFY). Lets start by eyeballing the plot of the stock prices of INFY for the past one year.

## Set the working directory using setwd() ##
# Reading the relevant file.
infy <- read.csv("01-10-2010-TO-01-10-2011INFYEQN.csv")

# Plotting the past one year's closing price of INFY
plot(as.Date(infy$Date, "%d-%b-%y"), infy$Close.Price, xlab= "Dates", ylab= "Adjusted closing price", type='l', col='red', main="Adjusted closing price of INFOSYS for past 1 year")

Eyeballing the above plot suggests that the series is NOT second order stationary. Meaning that the first two moments, of the distribution from which the data is drawn, changes with time. For a stationary series, the mean doesn't changes with time and the co-variance with any "k" lag is independent of "t" and it just a function of "k". But we see that both the conditions are violated above.

Let me attempt to explain the idea stationary in simple English language. For a moment suppose that you were to stand at time T = t and look at the value of the series, then look at the neighbors values to the left and right of "t", if by doing this exercise you can make out the value of "t" that you are standing at then it is possibly a non-stationary series. On the other hand if you were placed at time T = t in any stationary series, by doing the above exercise you would not be able to figure out the value of "t". (This definition came up during a discussion with Utkarsh some time ago).

A rule of thumb in any time series modelling is that we work with only stationary time series. If the series exhibits any non-stationarity, we have to remove that before we can employ any empirical analysis. In the above series the non-stationarity can be removed by using the returns instead of actual stock prices. (analogous to First differencing) .

## Calculating the returns of stock prices
infy_ret <- 100*diff(log(infy[,2]))

## Plotting the returns
plot(as.Date(infy$Date[-1], "%d-%b-%y"), infy_ret, xlab= "Dates", ylab= "Returns percentage(%)", type='l', col='red', main="Daily returns of INFOSYS for past 1 year")

We see that in the above plot the mean is fixed at 0 and the fluctuations are around that mean, that doesn't change with time. Now that we have taken care of the non-stationarity lets proceed on our task.

First we will plot the auto-correlation of the returns with the previous lags and see if there is any significant correlation that the returns have with the previous values.

## Plotting the ACF of INFY returns for the past one years
acf(infy_ret, main = "ACF of INFOSYS returns for past one year")

The blue dotted line is the 95% confidence interval. We can see that there is the 4th and the 7th lag significant in the ACF plot (there is one significant at 19th lag too but I choose to ignore that). Now lets see what I get if I regress the value of returns on the lagged values till lag 8th.

## Regressing the returns till the 7th lag
summary(lm(infy_ret[8:length(infy_ret)] ~ infy_ret[8:length(infy_ret) - 1] + infy_ret[8:length(infy_ret) - 2]+ infy_ret[8:length(infy_ret) - 3] + infy_ret[8:length(infy_ret) - 4] + infy_ret[8:length(infy_ret) - 5] + infy_ret[8:length(infy_ret) - 6] +infy_ret[8:length(infy_ret) - 7] ))## This is a simple OLS regression of the "inty_ret" starting from the 8th observation. I have started from the 8th observation to ensure that the number of obs. are same in the dependents and independent variables.

Output:
Coefficients:

                                 Estimate Std. Error t value Pr(>|t|)   
(Intercept)                      -0.09316    0.11321  -0.823  0.41140   
infy_ret[8:length(infy_ret) - 1]  0.08158    0.06479   1.259  0.20920   
infy_ret[8:length(infy_ret) - 2] -0.04017    0.06537  -0.614  0.53950   
infy_ret[8:length(infy_ret) - 3] -0.10049    0.06528  -1.539  0.12504   
infy_ret[8:length(infy_ret) - 4]  0.20153    0.06457   3.121  0.00203 **
infy_ret[8:length(infy_ret) - 5] -0.08566    0.06568  -1.304  0.19344   
infy_ret[8:length(infy_ret) - 6] -0.06849    0.06584  -1.040  0.29928   
infy_ret[8:length(infy_ret) - 7] -0.12395    0.06621  -1.872  0.06241 . 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Multiple R-squared: 0.08717, Adjusted R-squared: 0.05998

Only the coefficient of the 4th lag is statistically significant, and the Adjusted R-squared is a small 0.05998 (i.e ~ 6% of the explanation is provided by the above regression).

In the previous post we had reached the conclusion that the returns series is completely random (using runs.test()). But here we have fit in a model that provides ~ 6% of the explanation, the important question that needs to be addressed now is that the can we use this model to predict the stock returns (and make some money using a trading strategy that employs the above regression).

The model suggests that there is a statistically significant explanation that is being offered by the 4th lag in the above regression, but is this explanation economically significant? Now is when the economic intuition comes into play. The given sample data for the stock prices of INFY for the paste one year has confessed that the 4 days ago stock price provides a statistically significant explanation of today's stock prices. But a major point, perhaps the most important, that we are missing in the above model is the transaction costs or market micro-structures.

Meaning that a statistically significant 4th lag does not mean that the explanation offered is economically significant too. To check if the relation is economically significant, we will have to adjust the prices for transaction costs and then do the regression and see if we get a similar result. Efficient market hypothesis that this statistical significant will disappear once you account for these transaction costs (impact cost or cost of trading). It seems to be intuitive too, because if we look at the ACF plotted above the auto-correlations are not significantly different from 0 and once we account for the transaction costs the 95% band will also broaden.

So the lesson is that a simple regression of current returns on the lagged returns (auto regressive model in time series parlance) might not be a reliable trading strategy :-)

P.S. In case anyone wishes to replicate the exercise the data can be obtained from here.

Principal component analysis : Use extended to Financial economics : Part 2

2011-10-23T02:14:00.000+05:30

My previous post talked about how we can employ PCA on the data for multiple stock returns to reduce the number of variables in explaining the variance of the underlying data. But the idea was greeted with skepticism by many. A caveat to the application of PCA was that the meaning or the intuition behind the variable is lost in the computation of components, but since I was concerned with only explaining the variance in the underlying data it worked as a powerful method. This method was also adopted by Dr. Oyama in his research paper.

I came across this interesting research paper published in one of the most coveted journals, (International Journal on Finance and Economics) which has used the PCA on the macroeconomic variables. So relating it to my previous post and the question of investigating the factor affecting stock returns, what they do is that they take individual stocks and regresses it on the PCA of macroeconomic variables. At the first glance this would appear an inappropriate use of the tool, however the question that it intends to ask is also very different. There has been or infact still is this great battle between the Capital Asset Pricing (CAPM) school of thought and Arbitrage Pricing Theory (APT) school of thought. What the essential difference between these two schools is that the CAPM says that the returns to a stock are sensitive to only one factor, i.e the market rate of returns (single factor model) which captures the effect of all the various factors, however the APT guys scream that it depends on a number of factors and market rate of return would be just one of them. So in this battle to prove their point the APT folks take stock returns as a dependent variable and as independent variables they use the PCA of many macroeconomic variables. They do this essentially so that they can prove that there are more than one factors playing a significant part in the explaining the returns on a stock, they do not care whether the economic intuition is lost in the process on one-upmanship.

Well anything is fair in war, but I think this path taken by APT guys of proving their point is a fair point. The variables that are thrown into the PCA are chosen with an economic intuition in mind, so its not correct to say this entire methodology is flawed on the pretext that the PCA of these macroeconomic variables have no economics intuition. If one really digs into the intuition behind PCA as explained in this article one can visualize what the PCA of macroeconomics variables would be representing, the principal underlying components that give rise to such a macroeconomic series (signals). This is an unconventional way of doing an econometric study, as in this is more of a qualitative than a quantitative study. I cannot make sense of a statement that 1 unit increase in my first PCA results in "x" unit increase in returns, this is an absurd statement, but nevertheless the methodology of answering the underlying question is not absurd.

The methodology for calculating the PCA for macroeconomics is same as what was used in the previous post. We can have a number of macroeconomic variables, however, I have done this demonstration with only 2 variables change in Mumbai inter-bank offer rates (MIBOR) rates and change in INR/USD exchange rates.

#### Calculating the PCA of macroeconomic variables ####

# Reading the relevant files
mibor <- read.csv("MIBOR.csv", na.strings="#N/A")
exchange <- read.csv("Exchange_rates.csv", na.strings="#N/A")

# Making sure that there are no missing values in the data, the missing values are replaced by linear interpolation
mibor[, 2] <- approx(as.Date(mibor$Dates, '%d-%b-%y'), mibor[ ,2], as.Date(mibor$Dates, '%d-%b-%y'))$y # approx() returns the interpolated values in column 'y'

exchange[, 2] <- approx(as.Date(exchange$Year,'%d-%b-%y'), exchange[ ,2], as.Date(exchange$Year, '%d-%b-%y'))$y # Similarly for exchange rates

# Now we will have to compute the change in MIBOR and exchange rates:
for(k in 2:nrow(mibor))
{
mibor$Change1[k] <- ((mibor$MIBOR[k] - mibor$MIBOR[k-1])/mibor$MIBOR[k-1])*100
}

for(j in 2:nrow(exchange))
{
exchange$Change[j] <- ((exchange$Exchange.rates[j] - exchange$Exchange.rates[j-1])/exchange$Exchange.rates[j-1])*100
}

# Creating matrix of the data
macro <- as.data.frame(rep(0, 2498))

macro$ex <- exchange$Change
macro$rate <- mibor$Change1

macro.mat <- as.matrix(macro, nrow = 2498 , ncol = 2)

# Calculating the principal components

prin.macro <- princomp(macro.mat)

barplot(prin.macro$sdev[1:2]/prin.macro$sdev[1]) # bar plot of the ratio is the standard deviation of the PC's with the standard deviation of the 1st PC, this can help us decide which PC to be used by simple eyeballing the plot.

# Extract the 1st PC in a variable

load1 <- prin.macro$loadings[,1] # get the loadings of the 1st PC into 'load1'

pr.ma <- macro.mat %*% load1 # matrix multiply the loading with the original data to get the PC

pr1 <- as.numeric(pr.ma) # And 'pr1' has your first principal component.

Now if you have read this article that I shared earlier, the first component gives me the principal underlying factor that gives rise to the series (signals) that I observe.

Caveats and cautions:

Before I am bombarded with question about the relevance of this methodology in this particular problem let me confess that in the above particular exercise its not the best idea to use PCA. One should use PCA only if one has many macroeconomic variables like M1, M3, CPI, IIP, WPI, etc and there is high correlation between them. I have done the above only for illustrative purposes (since I am not employing PCA on macro-variables in my project).

Assuming I had other macro-variables in the calculation of my PCA above, I could use them as an independent variable and regress them on the stock specific returns to understand if there are more than one factors that come out to be significant in this regression. What this would essentially mean is that are there more than one factors explaining the stock returns. For the finance students, this is the classic argument between the Capital Asset Pricing (CAPM) guys and the Arbitrage Pricing Theory (APT) guys.

If in the regression of the PC's of macroeconomic variables on the returns we find that more than one factors are significant then we would have found evidence for or against one of the above theories.

Limitations in the above method:

In the above exercise I have just taken 2 variables merely for illustration, I would have to include more variables if I really intent to find evidences against the CAPM.
Again the caveat of interpretation applies. Its difficult to interpret the PCA on macrovariables

My take:

Despite of the above caveats I think PCA does a good job of giving us a good empirical way to answer the question of whether CAPM or APT is a superior asset pricing theory. The main idea behind these asset pricing theories is to see what factors need to be incorporated to compute the asset specific risk premiums, so on that front, its not difficult to see what people still stick to CAPM, as its a lot simple and easier to work with. APT on the other hand makes a bold assertion of multiple factors governing the returns but it has always found it difficult to list out these factors. Well in most cases you would be able to conclude that the returns can be modeled best using a multi-factor model, but whats interesting to ask is how much additional explanation is offered by factors other than the 1st one. If my multi factor model offers just a marginally higher adjusted-R-squared, is it reasonable for you to complicate your model is the question?

If it were only to prove this empirically why would you not choose the multi-factor model? The idea of modelling returns is to identify the factors that you need to take into account while computing the risk premiums, and if your factors do not have any economic intuition how can you assign risk premia looking at your principal component, that's the big question. But I would hate to impose my take here, further discussion is welcome.

P.S : In case you wish to replicate the above exercise you can find the data here, exchange rates, MIBOR.

Predictability of stock returns : Using runs.test()

2011-10-21T19:21:00.001+05:30

Financial market is interesting place, you find people taking positions (buying/selling) based on their expectations of what the security prices would be and are rewarded/penalized according to the accuracy of their expectations. The beauty of financial markets is that it provides a platform for everyone to come in with their respective expectations and allows them to interact and exchange securities. I emphasize on everyone because this everyone includes a auto-rickshaw driver, a clerk and also sophisticated econometricians and analysts. An obvious point then is that if your expectations are consistently correct, i.e you can predict the price movements before it happens on the exchange, you are a rich man. Assuming for all practical purposes that there is no oracle in our universe, who can do these predictions with 100% accuracy, the job of this prediction rests upon an econometrician/statistician. Lets see if they can do a good job too.

I took the stock returns data for INFOSYS (INFY on NSE) for the past one year and tried to see if I could make this data confess its underlying linear/non-linear generating process. I started by employing a rather simple, straight forward and easy to interpret Runs test. Its a non-parametric statistical test that will test the null hypothesis of whether the underlying series is identical and independent distributed. For those who are not too familiar with statistical parlance, non-parametric in simple term means that we have to make no assumptions about what the underlying data should be like. There is a huge surge in the applications of non-parametric statistics to explain various processes, this is because the biggest deterrence to conducting these kinds of tests, i.e the computational issues, are no longer a problem in this generation of rapid computation. The idea of empirical analysis is about trying to theorize a null hypothesis and then try your best to bring it down using empirical evidence. (analogous to Karl Popper's idea of falsification of a theory, you hang on to a theory so long as it has not betrayed you yet)

## Doing runs test on INFY daily returns
> infy <- read.csv("01-10-2010-TO-01-10-2011INFYEQN.csv") ## Reading the stock price data

> infy_ret <- 100*diff(log(infy[,2])) ## Since the second column in the data has the stock prices I have used [log(Pt) - log(Pt-1)]*100 as the returns.

> runs.test(factor(infy_ret > 0)) ## what this has done is that it has created a category variable that takes value 1 of infy_ret > 0 and 0 otherwise.

What this does is that tells me whether the runs of the returns are predictable, i.e say if I represent possitive return by + and negative return by - then my series of returns would probably look like +,+,-, +, -, -, -, +, ...
now that this test check is can I predict whether the next day will have + or -

Output:
Runs Test
data: factor(infy_ret > 0)

Standard Normal = 0.1308, p-value = 0.8959  ## High p-value means you cannot trash your null hypothesis.

For those not familiar with statistics, the p-value is nothing but the probability of you reject a null hypothesis when it is actually true. So in simple words it gives me the probability that I might end up rejecting a correct null hypothesis. (be very careful with the interpretation of p-value, many times people end up misunderstanding it, many a times even I have fallen prey to this). Therefore you cannot reject your null hypothesis under such a high probability of committing this error or wrongly rejecting a correct hypothesis , you just don't have enough evidence. Therefore your series is a random walk (you can understand this in the literal English language sense, but the definition is not so trivial in time series parlance).

P.S In case you want to replicate this exercise the data can be obtained from here.

Principal component analysis : Use extended to Financial economics : Part 1

2011-10-16T01:45:00.000+05:30

While working for my Financial economics project I came across this elegant tool called Principal component analysis (PCA)which is an extremely powerful tool when it comes to reducing the dimentionality of a data set comprising of highly correlated variables. This tool finds majority application in genetic research, which deals with data sets having many variables that are highly correlated.

I will try and be as explicit and refrain from using statistical/mathematical jargons to explain what/how about this tool . To state a few stylized facts PCA is used mainly for:

compressing the data
filter some of the noise in the data

Problem at hand:

I was trying to investigate the factors that affect the returns of stocks in the Indian equity market, however I wanted to take into account all the S&P CNX 500 companies. What would be really nice if I could somehow find a way of squeezing the 500 companies into say not more than 2-3 variables that can be representative of the entire set of 500 companies. This is precisely where PCA comes into play and does a fantastic job. What it gives me is just 1 variable that I can use instead of all the 500 companies!!!

Hats off and a bow of respect for the contributors/donors of packages to the CRAN servers that the above simplification can be achieved using just one line of script in R. Sounds easy, but what one really needs to do is to understand what PCA does and how the output from this script can be interpreted. Again, at the risk of over simplication (however trying hard to maintain my commandment of simplicity), I would illustrate in a crude manner the working of PCA.

What PCA does:

Let me explain this relating to the above example, if I do a PCA on the returns data for the 500 companies, I would obtain 500 principal components. These components are nothing but the linear combination of the existing 500 variables(companies) arranged in the decreasing order of their variance. So 1st principal component (PC)has the maximum variance and 500th principal component (PC)has the least variance. The variance in the PCA represent nothing but the variance in the data. So 1st PC explains the maximum amount of variance in my data. One magical feature of PCA is that all these 500 components will be orthogonal to each other, meaning these components will be uncorrelated with each other. So essentially if we look at PCA as a black box it takes inputs as data set of highly correlated variables and gives as output PC's that explain the variance in the input data and they are uncorrelated with each other.(I don't leverage this feature in this particular problem, I would illustrate this use in other part of this blog)

How PCA does it:

Since I have taken a vow of simplicity, I dont have much to say here.:-) However for the mathematically inclined and certainty freaks like Madhav, this paper does a brilliant job of illustrating the matrix algebra that goes behind PCA computations. There are essentially 2 methods of calculating PCA, one is the eigenvalue decomposition (done using princomp() command in R)and the other is singular value decomposition (done using prcomp() command using R).

How this can be done in R:

####### Calculating Principal component of returns of S&P CNX 500 companies ########

## Access the relevant file ##

returns <- read.csv("Returns_CNX_500.csv")

One caveat that you need to keep in mind in that there should be no "NA" values in your data set. A presence of an NA would impede the computation of the var-covar matrix and hence their eigen vectors(i.e the factor loadings)

## Dealing with missing values in the returns data for companies

for(i in 2:ncol(returns))

{

returns1[, i] <- approx(returns$Year, returns1[ ,i], returns$Year)$y ## approx function basically fits the value of linear approximate between the missing data points and the column $y stores the approximated values.

}

## Convert the data into matrix ##

ret <- as.matrix(returns1, nrow = dim(returns1)[1], ncol = dim(returns1)[2])

##Computing the principal component using eigenvalue decomposition ##

princ.return <- princomp(ret) ## This is it.!!

## Identifying what components to be used ##

barplot(height=princ.return$sdev[1:10]/princ.return$sdev[1]) ## I am plotting the standard deviation of the PC's divided by standard deviation of PC 1, this can help us decide on a benchmark that we can use to select the relevant components.

Standard deviation of the first 10 components compared to 1st PC

We can clearly see from the above figure that as expected the first PC does the majority of the variance explanation in the returns data for the 500 companies. So if we want to identify factors that influence the returns of S&P CNX 500 companies I can use the 1st PC as a variable in my regression. So far we have calculated the principal components, now we will extract out 1st PC as a numeric variable from the matrix.(princ.return)

## To get the first principal component in a variable ##

load <- loadings(princ.return)[,1] ## loadings() gives the linear combination by which our input variables will be linearly weighted to compute the components, and this command gives us the loading for 1st PC.

pr.cp <- ret %*% load ## Matrix multiplication of the input data with the loading for the 1st PC gives us the 1st PC in matrix form.

pr <- as.numeric(pr.cp) ## Gives us the 1st PC in numeric form in pr.

One question that might be raised is why not just use the S&P CNX 500 index returns as an input to the regression? The simple answer to that question would be that PC 1 gives you a relatively clear signal of the returns as opposed to the index which would have a lot of noise. This question would have made sense in the 1900's when the technology was not so efficient in terms of computation. Since now computational time and effort finds minimum weight in any researchers mind there is no reason to settle for anything but the best.

There is an important caveat that must be kept in mind while doing analysis using PCA, though PCA has a clear mathematical intuition it lacks an economic intuition. That is, one unit change in PC 1 of returns has a mathematical meaning but no economic meaning, you cannot make sense of this statement that PC 1 of returns for the 500 companies has gone up by "x" amount. Therefore the use of this analysis should be limited to factor analysis and not to be extended to predictive analysis.

In case you wish to replicate the above exercise the data can be obtained from here.

Modelling with R: part 4

2011-10-13T15:38:00.002+05:30

In part 3, we ran a logistic model to determine the probability of default of a customer. We also saw the outputs and tried to judge the performance of the model b plotting the ROC curve. Let's try a different approach today. How about a decision tree? For starters, we being a simple decision tree with no priors.

##########-------------2.2: Recursive partitioning------------###############

## 2.2.1: Decision tree with no priors ##

What is recursive? How is partitioning done? How does partitioning translate to predictions? What are priors? Let's try and decode these big words that we have written above.

Partitioning means splitting the data into subsets based on attributes. Partitioning typically happens in a top-down approach by choosing a variable, or an attribute, that is the "best" attribute to split the data. This "best" variable would be one that splits the data into homogeneous subsets with each subset having the same category of the target variable, i.e, either all zeros or all ones. It is in the way of choosing this "best" variable that different tree construction algorithms differ. For example, CART, uses Gini coefficient to choose the best variable. I read that ID3 and C4.5 use the concept of entropy, but I really don't know how these function.

After a subset is formed, the splitting process continues in a "recursive" manner and the recursion halts when the subset at a node has the same value for the target variable.

To run a decision tree we need to call the "rpart" library.

library(rpart)

We use a recursive partioning tree model to predict the likelihood of default.

t1.noprior <- rpart(default ~

amt.fac +

age.fac +

duration +

chk_acct +

history +

purpose +

sav_acct +

employment +

install_rate +

pstatus +

other_debtor +

time_resid +

property +

other_install +

housing +

other_credits +

job +

num_depend +

telephone +

foreign

, data = dev)

plot(t1.noprior)

# Plots the trees

text(t1.noprior)

# Adds the labels to the trees.

Before we score the validation data, we need to prune the tree to ensure that we haven't over fitted the data.

In order to prune the data we need to find the appropriate complexity parameter. The ideal value of the complexity parameter in this case would then be the one that minimizes the X-error.

We can find this value as

printcp(t1.noprior)

# Here we see that the value is 0.013514.

We can write a small script to find this value

t1.noprior.cp <- t1.noprior$cptable[which.min(t1.noprior$cptable[, "xerror"]), "CP"]

t1.noprior.prune <- prune(t1.noprior,t1.noprior.cp)

We need to score the pruned tree model the same way we did for the Logistic model.

val$t1.yhat <- predict(t1.noprior.prune, val, type = "prob")

Now notice the difference between the two predictions- one of the Logistic model and the other of the tree model. The logistic model predictions have only column, which gives the probability that the person will default, or mathematically, the probability of y being equal to 1. On the other hand the tree model predictions give two columns, first for y being equal to zero and second column for y being equal to 1. Hence, for plotting the ROC curve we need to take only the probability of y being equal to 1.

We will also plot the ROC curve for this tree.

t1.scores <- prediction(val$t1.yhat[,2], val$default)

t1.perf <- performance(t1.scores, "tpr", "fpr")

# Plot the ROC curve

plot(t1.perf, col = "green", lwd = 1.5)

Now how do we compare the performance of this decision tree model with our earlier logistic model. One simple way would be to plot the two ROCs in the same graph and visually guess it.

# Add the ROC curve of the logistic model and the diagonal line

plot(m1.perf, col = "red", lwd = 1, add = TRUE)

abline(0, 1, lty = 8, col = "grey")

legend("bottomright", legend = c("tree", "logit"), col = c("green", "red"), lwd = c(1.5,1))

Similarly we can calculate the KS statistic as well

ks1.tree <- max(attr(t1.perf, "y.values")[[1]] - (attr(t1.perf, "x.values")[[1]]))

ks1.tree

Another way to compare the performance of two models, especially when the dependent variable is binary, is the Area Under the Curve or the "AUC"

# AUC for the logistic model

m1.auc <- performance(m1.scores, "auc")

m1.auc

# AUC for the decision tree

t1.auc <- performance(t1.scores, "auc")

t1.auc

Modelling with R: part 3

2011-10-05T18:40:00.001+05:30

The previous posts, part 1 and part 2, detailed the procedure to successfully import the data and transform the data so that we can extract some useful information from them. Now it's time to get our hands dirty with some predictive modelling.

The dependent variable here is a binary variable taking values "0" and "1", indicating whether the customer defaulted or not. There are specific, but many, modelling techniques that help us with binary dependent variables. Let's start with the one of the simplest, the Logistic Regression model. And this time, let's just dive straight into it and try to make sense as we go along.

##########-------------2.1: Logistic regression------------###############

m1.logit <- glm(default ~
amt.fac +
age.fac +
duration +
chk_acct +
history +
purpose +
sav_acct +
employment +
install_rate +
pstatus +
other_debtor +
time_resid +
property +
other_install +
housing +
other_credits +
job +
num_depend +
telephone +
foreign
, family = binomial(link = "logit"), data = dev)
# The glm command is for "Generalized Linear Models"
# "~" is the separator between the dependent (or target variable) and the independent variables
# Independent variables are separated by the "+" sign
# "data" requires the data frame in the which the variables are stored
# "family" is used to specify the assumed distribution.

The above code is pretty much self-explanatory. We are running a logistic model on our transformed data. We have specified the "logit" link which means that the estimates (or coefficient or weights) will be for the log of odds ratio and not the direct probabilities. Additionally, note that in case we had to regress the dependent variable, "default" on all the independent variables, we could have simply written
# m1.logit <- glm(default ~ ., family = binomial(link = "logit"), data = dev)
Here the "." indicates all the variables in the data set except for the dependent variable. We are not using the command here is because we have the original "amount" and "age" variables as well as the transformed "amt.fac" and "age.fac" variables as well and we should use only one form, either the continuous form or the factor form but not both.

We can check the output of the model as
summary(m1.logit)

You will see an output that looks like this
Call:

glm(formula = default ~ amt.fac + age.fac + duration + chk_acct + 
    history + purpose + sav_acct + employment + install_rate + 
    pstatus + other_debtor + time_resid + property + other_install + 
    housing + other_credits + job + num_depend + telephone + 
    foreign, family = binomial(link = "logit"), data = dev)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.9970  -0.7118  -0.3821   0.6864   2.5393  

Coefficients:
                  Estimate Std. Error z value Pr(>|z|)    
(Intercept)        0.36717    1.28298   0.286  0.77474    
amt.fac2600-5000  -0.39145    0.28742  -1.362  0.17322    
amt.fac5000+       0.25237    0.40440   0.624  0.53258    
age.fac30-40      -0.12406    0.26915  -0.461  0.64486    
age.fac40+        -0.40640    0.30975  -1.312  0.18950    
duration           0.03454    0.01152   2.998  0.00272 ** 
chk_acctA12       -0.26207    0.26204  -1.000  0.31725    
chk_acctA13       -0.71081    0.43527  -1.633  0.10246    
chk_acctA14       -1.66798    0.27841  -5.991 2.08e-09 ***
historyA31         0.51820    0.68122   0.761  0.44684    
historyA32        -0.44583    0.54524  -0.818  0.41354    
historyA33        -0.36003    0.58949  -0.611  0.54137    
historyA34        -0.89606    0.54235  -1.652  0.09850 .  
purposeA41        -1.30443    0.44483  -2.932  0.00336 ** 
purposeA410       -0.68765    0.79700  -0.863  0.38825    
purposeA42        -0.50703    0.31109  -1.630  0.10314    
purposeA43        -0.95899    0.29629  -3.237  0.00121 ** 
purposeA44        -1.07335    0.95318  -1.126  0.26013    
purposeA45        -0.40166    0.64852  -0.619  0.53569    
purposeA46         0.07051    0.48321   0.146  0.88399    
purposeA48        -0.79836    1.32004  -0.605  0.54531    
purposeA49        -0.63881    0.40523  -1.576  0.11493    
sav_acctA62       -0.20567    0.35346  -0.582  0.56064    
sav_acctA63       -0.29441    0.47750  -0.617  0.53752    
sav_acctA64       -2.03887    0.74864  -2.723  0.00646 ** 
sav_acctA65       -0.95937    0.32147  -2.984  0.00284 ** 
employmentA72      0.24498    0.51552   0.475  0.63463    
employmentA73      0.01948    0.50334   0.039  0.96913    
employmentA74     -0.41741    0.54075  -0.772  0.44016    
employmentA75      0.01317    0.50928   0.026  0.97936    
install_rate       0.22383    0.10671   2.098  0.03594 *  
pstatusA92        -0.71859    0.53450  -1.344  0.17882    
pstatusA93        -1.19509    0.52355  -2.283  0.02245 *  
pstatusA94        -0.95261    0.60277  -1.580  0.11402    
other_debtorA102   0.71012    0.48442   1.466  0.14267    
other_debtorA103  -0.84087    0.51928  -1.619  0.10538    
time_resid         0.06092    0.10684   0.570  0.56851    
propertyA122       0.43538    0.31755   1.371  0.17035    
propertyA123       0.45504    0.29062   1.566  0.11740    
propertyA124       0.95268    0.52675   1.809  0.07051 .  
other_installA142  0.54490    0.47279   1.153  0.24911    
other_installA143 -0.45903    0.28664  -1.601  0.10929    
housingA152       -0.58427    0.29480  -1.982  0.04749 *  
housingA153       -0.95369    0.59035  -1.615  0.10621    
other_credits     -0.08412    0.24145  -0.348  0.72753    
jobA172            0.51201    0.80390   0.637  0.52419    
jobA173            0.56316    0.77175   0.730  0.46556    
jobA174            0.09850    0.77019   0.128  0.89823    
num_depend         0.38802    0.29445   1.318  0.18758    
telephoneA192     -0.19916    0.24766  -0.804  0.42131    
foreignA202       -1.35004    0.85752  -1.574  0.11541    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 853.51  on 699  degrees of freedom
Residual deviance: 628.93  on 649  degrees of freedom
AIC: 730.93

Number of Fisher Scoring iterations: 5

We can see that we have five variables which are significant at 1 percent level of significance and one variable significant at 0.1 percent level of significance. Great! Or not so great! How do we judge how good/bad our model is? Have no fear, R is here.

One of the most popular methods to check the performance of a model, which has a binary independent variable, is by plotting the Receiver Operating Characteristic (ROC) curve which is a plot of the sensitivity, or the true positive rate vs. the false positive rate. In simpler terms, it plots the no. of 1s in the model correctly identified as 1 vs. the no. of 0s in the model which are incorrectly identified as 1.

Also, it makes more sense to plot the ROC curve on the validation sample. For this, we would need to score that validation data set. By "scoring" we mean that with the values of the independent variables given in the validation sample, we will run the model equation on each row of the data set and get the probability of the customer defaulting. Simply, use all the Xs and the Y-hats or the predicted Ys.
val$m1.yhat <- predict(m1.logit, val, type = "response")
# The predict command runs the regression model on the "val" dataset and stores the estimated y-values, i.e, the yhat.

library(ROCR)
m1.scores <- prediction(val$m1.yhat, val$default)
# The prediction function of the ROCR library basically creates a structure to validate our predictions, "val$yhat" with respect to the actual y-values "val$default"

We can then plot the ROC curve by plotting the performance fuction which checks how our model is performing.
plot(performance(m1.scores, "tpr", "fpr"), col = "red")
abline(0,1, lty = 8, col = "grey")
# "tpr" and "fpr" are arguments of the "performance" function indicating that the plot is between the true positive rate and the false positive rate.
# "col" is an argument to the plot function and indicates the colour of the line
# "abline" plots the diagonal, "lty" is the line type which is used to create the dashed line

The ROC curve's performance is judged by the area under the curve formed by the ROC curve and the 45' line. The idea is that if we randomly try to classify the customers who will default, we will get the 45' line. And if we use a predictive approach, then the approach should try to classify all the 1s and 1, i.e, correctly identify all the people who defaulted. Such a curve would then be a right triangle with vertices at (0,0), (0,1), (1,1). But since predictive modelling is usually accompanied with some error, we get a curve like we see in the graph. And the rule of thumb is, more the area under the curve, better the model.

In addition to the ROC curve, the KS statistic is also a popular metric used to judge model performance, the Kolmogorov-Smirnov statistic. It is the maximum difference between the cumulative true positive rate and the cumulative false positive rate

For this let's store the output of the "performance" function above to an R object
m1.perf <- performance(m1.scores, "tpr", "fpr")
ks1.logit <- max(attr(m1.perf, "y.values")[[1]] - (attr(m1.perf, "x.values")[[1]]))
ks1.logit

Well, the ROC shows that the model is not so great but not so bad either. So what do we do to improve our model? Do some more transformations, tweak it a bit, or just try a different technique?

I'll leave this as an open question to make it interactive. Let's try all the "feasible" approaches that people suggest that will help us improve the model.

Simple time series plot using R : Part 2

2011-10-04T22:30:00.002+05:30

I would like to share my experience of plotting different time series in the same plot for comparison. As an assignment I had to plot the time series of Infant mortality rate(IMR) along with the SOX emission(sulphur emission) for the past 5 decades in the same graph and compare how the intensities have been varying in the past 5 decades.Well to start with there is a problem of how to get these plots in the same graph as one is a mortality rate and other is an emission rate, having different units of measurements!

What we essentially want to do is to see how the intensity of these problem has been changing, whether the intensity has increased/decreased. From a policy makers standpoint, whether it requires immediate attention or not. So what we can do instead is divide all the IMR values by the maximum IMR value that we have for the past 5 decade and store them as "IMR.std". Similarly divide all the SOX values by the maximum SOX value and store it in "SOX.Std". What we have achieved is a parsimonious way of representing the 2 variables that have values between "0 and 1".(Achieving the desired standardization).

Now that we have "IMR.Std." and "SOX.Std." both with values between "0-1" I can plot them in the same graph. Recalling from the previous post :

# Make sure the working directory is set to where the file is, in this case "Environment.csv":

a <- read.csv("Environment.csv")

# Plotting the "IMR.Std" on a graph

plot(a$Year, a$IMR.Std., type="l", xlab= "Years", ylab= "Intensity of problem(normalized between 0-1)", col="green" , lwd=2)

# Adding the plot for "SOX.Std."
# lines(...) command basically adds the argument variable to the existing plot.

lines(a$Year, a$SOX.std., type="l", col="black", lwd=2)

Ideally this should have done the job. Giving me the IMR.Std(in green) and SOX.Std.(in red) in the same plot. But this dint happen for the reason that the data for SOX was available only after 1975 and also the data was not available for alternate years. Well I thought that R would treat it trivially and just plot the "non-NA" values of SOX.std. that were there but as it happens that this was not such a trivial thing for R. It demands a lot more rigor(just like a mathematical proof), to execute a command, not taking anything for granted. Hence to get the desired result I had to specify that it considers only the "non-NA" values for SOX.Std.

The code for SOX.std had to be altered a bit:

a <- read.csv("Environment.csv")

plot(a$Year, a$IMR.Std., type="l", xlab= "Years", ylab= "Intensity of problem(normalized between 0-1)", col="green" , lwd=2)

# All I need to make sure now is that I direct R to refer to only the "non-NA" values in the SOX.Std. variable.
lines( a$Year[ !is.na(a$SOX.std.) ], a$SOX.std.[ !is.na(a$SOX.std.) ], type="l", col="black", lwd=2)

# If there are 2 separate graphs on the same plot to label the different graphs(depending on their color in this case) we use the legend(...) command :

b <- c("IMR rate", "SOX")
legend("bottomleft", b , cex=0.8, col=c("green", "black"), lwd=2, bty="n")

# title(...) command gets a main title to the plot
title(main="Stage model for country X", col.main="Black", font.main=4)

This is how the plot finally looks like.

It was Utkarsh's generosity, who gave me the codes, that saved me a lot of time in solving this small issue, I wish to pass this on as it might save someone else's.

Modelling with R: part 2

2011-10-02T19:39:00.002+05:30

I apologize for the delay in the second post (just in case anybody was waiting), I had been vary involved with work the past week. I shall try to be more regular. Well, in the previous post, we successfully imported data into R and got a basic "feel" of it by looking at the various variables present and their types as well. Now we will try to process the data to make sense of it. Data by themselves are just space hogging particles, aesthetically challenged, and practically worthless... unless... well... unless... we can get some information out of them. And to get that information, they need to be processed, transformed, and at times coerced. This post will describe how we can start to do that. So, let's grab its throat and make it spit out the ugly truth (excuse me for the histrionics). Also, this last sentence bore no relation this horrible Gerad Butler movie.

##########----------------1.2: Processing the data------------###############

There are three main steps involved here:

1. Preliminary visulaization

2. Data transformation and/or variable creation

3. Development-validation division of data set

Let's start with 1.

Suppose we want to check the distribution of amount in the given data. We can use a simple plot command.

plot(amount, type = "l", col = "royalblue")

plot(age, type = "l", col = "brown")

Now the plots we have created present the variation in these variables in a manner which neither easily discernible nor is it clean.

It would be better if created histograms to check the frequency distribution of these variables.

hist(amount)

hist(age)

# The "hist" command has a lot of options that help extend the features of the plot.

Now, with the hist command, we are able to see the picture clearly (literally). But it is not always appropriate to plot a histogram. What is the best way depends largely on the problem at hand. Suppose we had a multiple time series and we wanted to check the behaviour, then it would be better to use the plot command which will not only present the data in a much neater way but also enable us to compare different time series in a single plot. For a simple example, you can check Shreyes' post.

Coming back, we can similarly visualize the pattern, frequency and distributions of other variables as well.

Now, in case you have ever worked on a credit scoring exercise before, you might have heard that it is better to create categories out continuous variables. This helps a lot while implementing the model that we build because it is more convenient to come with strategies for individuals belonging to a particular income group rather than for all individuals with specific incomes.

For this we need to bin some variables like amount and age. One approach to do this is to run the following code

DO NOT run this chunk of code. I will explain later why.

# g.data$amount <- as.factor(ifelse(amount <= 2500, "0-2500",

ifelse(amount <= 5000, "2600-5000", "5000+")))

Here we are creating three categories for the variable "amount". One for those with income level less than 2500, one for income between 2500 and 5000, and one for income greater than 5000.

There is an important point to note here. Above, while creating the category variable, we overwrote the original variable "amount" in the R object "g.data". Ideally this process is not well advised because if we later find that there was an error in our code or there was some flaw in the logic and we need to change it, we will have to re-do all the steps that we have done till this point. But, there is another side here as well. R, while working, stores all the data and the objects that we create in the RAM and hence if the data set is of considerably large size then creating additional variables by transformation is not a very wise idea either. This trade-off needs to be balanced.

In this case, the data set is quite small and hence it would be better if we create an additional

object instead of overwriting the original one.

g.data$amt.fac <- as.factor(ifelse(amount <= 2500, "0-2500",

ifelse(amount <= 5000, "2600-5000", "5000+")))

head(g.data$amt.fac)

Similarly, we can do so for "age".

g.data$age.fac <- as.factor(ifelse(age<=30, '0-30', ifelse(age <= 40, '30-40', '40+')))

head(g.data$age.fac)

Here our dependent variable is "response". It is a factor variable and has "1" and "2" as the factor levels. Now, R by itself can handle factor variables and so we do not need to transform them unless we plan to combine categories. But I like to keep the response category coded as "1" (this is just because of habit and nothing else). Hence, I reassign the levels to "0" and "1".

g.data$default <- as.factor(ifelse(response == 1, "0", "1"))

is.factor(g.data$default)

contrasts(g.data$default)

head(g.data$default)

We attach the data again to include the newly created variables.

attach(g.data)

In the previous post, one the comments introduced me to the with() command as a substitute for the attach(). The with() command also serves the purpose quite well. It reduces the pain of writing the object name with the $ sign before we can refer to a variable but it needs to be included for every operation that we perform on the object.

Now, we saw that there are a lot of categorical variables present in the data. R provides many functions to plot categorical data.

Let's see an example.

mosaicplot(default ~ age.fac, col = T)

mosaicplot(default ~ job, col = T)

mosaicplot(default ~ chk_acct, col = T)

We can also use a spine plot.

spineplot(default ~ age.fac)

We can also check the relations between variables

library(lattice)

xyplot(amount ~ age)

In case you don't have the "lattice" library installed, you can download it by running

install.packages("lattice")

We can also condition on a variable and see the interaction

xyplot(amount ~ age | default)

"lattice" package also has the option for a barchart and it lets you plot the barchart and a histogram for factor variable type as well.

barchart(age.fac, col = "grey")

barchart(amt.fac, col = "grey")

histogram(employment, col = "grey")

histogram(sav_acct, col = "grey")

As a last step in this stage, we need to create a development sample and a validation sample. We take about 70% percent of the data as development sample and 30% as validation sample.

d <- sort(sample(nrow(g.data), nrow(g.data)*0.7))

# The sample command here creates a random sample of the number of rows in "g.data" and then

# selects 70% of this random sample and assigns it to object "d".

Note that here the sample is being generated from the row numbers and not the exact rows of data

so that if you see the object "d", you will see 700 randomly selected natural numbers between 1 and 1000 which are nothing but the row numbers in the data frame "g.data".

# The "sort" command in the beginning just sorts these randomly generated row numbers in an ascending order.

Then to create the development sample, we use the vector properties of R and assign the "d" rows to the R object "dev", and the remaining to the R object "val".

dev<-g.data[d,]

val<-g.data[-d,]

After creating the sample, we can check the size of the two samples vis-a-vis the original data.

dim(g.data)

dim(dev)

dim(val)

Finally we have been able to domesticate the data. We have sliced and diced them according to our needs. In the next post we will try to cook them in the modelling pan.

Modelling with R: part 1

2011-09-30T03:32:00.000+05:30

When I started work about 3 months ago, I didn't know much more than loading data and executing standard Econometric commands in R. But now I feel much much much more confident in using R for work, for research, for puzzles, and sometimes just for fun. I have learnt a lot of R, statistics, and real world application of Econometrics than I ever did at school (this in spite of the fact that I graduated from a school which has a very strong foundation in the field).

Much R that I have learnt has come from various R enthusiasts who have shared their knowledge and creativity of this elegant language over the internet. And as I grow with R, I found it quite compelling to share my learning journey as well. So, I have decided to describe some of the modelling techniques that are available in R. Along with this, I shall keep giving a few insights regarding the foundation and the need for the technique applied. Since all the techniques are quite big for one blog post, I have decided to divide entire journey into small, channelized, understandable and presumably fun exercises that will help the reader navigate through the models easily and efficiently.

For illustrating these techniques, I have taken the quite popular German credit data set which can be downloaded from here and here. Additionally, most of the techniques that I describe here are taken from this excellent guide. These blog posts can be thought of an extension to this guide. The German credit data has 1000 rows and 21 columns including the dependent variable, which in this case is binary- 1 means "good credit" and 2 means "bad credit". We need to predict whether a given case example will be a "good credit" or a "bad credit".

Now, there are a lot of techniques that can be used for such cases which have a dichotomous dependent variable. We will go through some of them in the following posts. For starters, here we will just try to import the data successfully and try to understand what kind of variables do we have.

###########================================================###############

# Section 1: Data preparation #

###########================================================###############

##########----------------1.1: Reading the data------------###############

Set the directory to where the german data file is located

setwd("D:/Softwares/R/Training/")

Read the data by importing the csv file

g.data <- read.csv("german.data.csv", header = F)

Note that by default the read.csv file has header the option set to TRUE. In this case however, since we don't have the variable names and hence the header row, we need to change the default option and set it to FALSE otherwise the first row of data will be taken as the header row.

# There is a slight mistake in the above code. I came to realize this when it was pointed out in one of the comments. As it happens, the data set available at the two links above is in a "space" delimited format. So the usual csv command will not work. I imported it as a csv because I had converted the tab delimited format to csv using a spread sheet (quite "uncool" I know, but I didn't know much R then). Anyway, thanks to another comment, we can use this command to import the file.
g.data <- read.delim("german.data",header=F,sep=" ")

It's always good practice to check the first few observations and see that the data were read in correctly.

head(g.data)

# "head" displays the first six observations

We can and should also check the number of rows of data and the number of variables present.

dim(g.data)

Now, in order to make our analysis easier to understand and implement, it is advisable to identify all the variables by specifying the column names. The column names can be found in the data dictionary.

names(g.data) <- c("chk_acct", "duration", "history", "purpose", "amount", "sav_acct", "employment", "install_rate", "pstatus", "other_debtor", "time_resid", "property", "age", "other_install", "housing", "other_credits", "job", "num_depend", "telephone", "foreign", "response")

head(g.data)

# The "names" function gives the column names of the data frame. Here we are assigning names using this function in a serial manner.

Well, if I may use some corporate jargon, it is advisable to do a DQ and DI check. For starters, it is important to find out the characteristic of all the variables present in the data. We can start with checking the type for each variable. For variables that we expect to be numeric or factor we can check them as

is.numeric(g.data$property)

is.factor(g.data$property)

is.numeric(g.data$age)

is.double(g.data$amount)

is.numeric(g.data$amount)

This process becomes quite a tedious exercise if we have a large number of variables R can make the task easier for us.

str(g.data)

# "str" will compactly display the structure of any R object

Additionally, to make the code more legible and easier to write, we can attach the object "g.data".

Then we can refer to g.data$amount simply as amount. But we need to ensure that there are not more than one objects with same name, i.e, if you have say two data sets data1 and data2, and both the data sets have the variable "amount", then if we attach both the data sets and call "amount", then "amount" will refer to the former object we attached.

attach(g.data)

head(amount)

Well, this completes the data import. We have successfully manged to import the data and correctly identify all the variables present. We will continue playing with data in the next post.

Project Euler: problem 6

2011-09-27T20:07:00.001+05:30

The sum of the squares of the first ten natural numbers is,

1² + 2² + ... + 10² = 385

The square of the sum of the first ten natural numbers is,

(1 + 2 + ... + 10)² = 55² = 3025

Hence the difference between the sum of the squares of the first ten natural numbers and the square of the sum is 3025 - 385 = 2640.

Find the difference between the sum of the squares of the first one hundred natural numbers and the square of the sum.

This is one quite simple.

# Create a vector with the first hundred natural numbers

x.1 <- (1:100)
# Square each element of the vector of natural numbers, i.e, square each natural number and store in another vector
y.1 <- x.1^2
head(y.1)
# Sum the first hundred natural numbers and the squares of each of the first hundred natural numbers and take the difference of these sums
a.1 <- sum(y.1)
b.1 <- (sum(x.1))^2
b.1 - a.1
b.1

Answer: 25164150

After solving a couple of these problems, andafter reading some solutions posted by aatrujillo here and here, I realize that my solutions are notgeneral, i.e., they only cater to the problem at hand and hence their scope isvery limited. For example, consider the above problem. Had the question askedto do the same analysis on the first 200 natural numbers, I would have torewrite the entire loop again. I understand that in this case it does notinvolve much more than changing the size of the x.1 vector, but for a problemthat involves more than one loop, it seems to be very "uncool". As aresult, I have decided to orient my results towards general solutions and thensolve the problem by specifying the parameters. Let's see how that goes. :)

Simple time series plot using R : Part 1

2011-09-21T23:15:00.001+05:30

As a task for my Financial eco assignment I had to plot a simple time series of the overnight MIBOR(Mumbai interbank offer rates) for the past one year . The job could very well have been done easily in MS-Excel but I choose to plot it in R instead and the quality of the graph, pixel-wise and neatness wise, was way better than what I could have obtained with MS-Excel. All this at the cost of a minimal 3 lines of code:

# The overnight MIBOR rates were stored in a file name "Call_Rates_2011.csv", this is just a normal Excel file saved in a CSV(comma separated delimited) format that R can read.
# The way in which R conceptualizes the data is similar to that in Excel, to draw a simple analogy you can assume that the variable "a" now stores the entire Excel spreed sheet in it.
# You will have to make sure that the working directory is the one that contains the file "Call_Rates_2011.csv"

a <- read.csv("Call_Rates_2011.csv")

# The 2 column headers in my CSV file were "date" and "mibor", so the below code plots "date" on the x-axis and "mibor" on the y-axis. The as.Date() tells R that the column "date" contains dates in the format "day-month-year"('%d-%b-%y').
# a$(column header) is the standard way of referring to a column in the "spreadsheet" contained in "a"
# xlab : x-axis label
# ylab : y- axis label
# type : line(l)
# col : color of the line

plot(as.Date(a$date,'%d-%b-%y'), a$mibor, xlab= "Months", ylab= "MIBOR overnight rates(percentage)", type='l', col='red')

# This is to get the titles in place
# main : main title
# col.main : color of the main title
#font.main : font size of the title

title(main="Overnight MIBOR rates for last one year", col.main="black", font.main=4)

And the plot hence obtained thus looks like:

Incase you can't make out the difference in the quality of the plot obtained just drop in your comments and email address and I will mail you the pdf and the jpg image of the plot. You can pull/stretch it to see that the pixels don't get distorted and it looks way neater if you present it in your slide in a presentation.

Project Euler: problem 3

2011-09-21T22:01:00.001+05:30

The prime factors of 13195 are 5, 7, 13 and 29.
What is the largest prime factor of the number 600851475143 ?

This one was quite easy, and much easier in R as it turns out.

The GNU Multi-Precision Library (GMP) is available as a package in R. So the only thing I had to do is install the library. Rest... well... not much...

library(gmp)
factorize(600851475143)

# The factorize function lists down all the prime factors of the number in the parentheses in ascending order.

Answer: 6857

Well, I understand I did not do much here than just being aware of something called the GMP. Don't blame me, blame them or him. But in due respect for the only concrete language known to humankind (Math not R), I shall try to come up with a more genuine approach.

Introduction to Beamer

2011-09-17T18:19:00.001+05:30

A friend of mine, who is quite smart by the way (she is a PhD. student in Physics at Cambridge), recently asked me for some help with Beamer. Well most of my knowledge and code came from Utkarsh when I had started about a year ago. Initially, I had just asked him for some help with Ubuntu, then a little bit with R, and then some with LaTeX, followed a bit by Python. This was all accompanied with constant questions from my side regarding mathematics, programming, web scarping, regular expressions etc. Most of the times he ended up writing the entire code for me. Though this was more or less like spoon-feeding a kid, he made it a point to ensure that I understood what the code was doing. After writing a lot of code for me and fixing a lot of bugs in the code that I wrote, I finally started to get some confidence while writing my own scripts. A year later I realize that much of teaching much did not emphasize that I get the "grammar" of script right, but that I got the "approach" of the script right. To be more clear, it's more important to know exactly what we want are trying achieve through the code rather than the way to do it. If you know what you want, you'll figure out a way to get it.

Well, when my friend asked for some help with Beamer, I didn't know how to exactly go about it. Should I forward her some documentation that I have stored in my computer which I never read myself? Well, I did do that. But the idea was a bit discomforting. Also, all the time that Utkarsh spent with an idiot like me should not go waste. So I tried to recollect some of the things that he taught me, scrapped through some of my beamer presentations that I had made in college and organized it into a TeX document.

Hope it helps...

% This is a very humble attempt to introduce fellow students to amazing world of LaTeX, and more specifically Beamer.

% This by no means is a guide or a reference manual for Beamer. It is simply a beginning. If you are planning your first presentation, this should serve the purpose and at the same time arouse curiosity to delve deeper.

% Also, I would like to point out that most of the code here has been acquired from various secondary sources and are not

% my creation. I have just assembled some information that I have into a organized document. I sincerely thank all those creative individuals you were generous enough to share their solutions with the world.

\documentclass[]{beamer}

\usepackage{graphicx}

\usepackage{amsmath}

\usepackage{tabls}

\usepackage{booktabs}

\usepackage{float}

\usepackage{array}

\usepackage[english]{babel}

\usepackage{hyperref}

\usepackage{subfigure} % This package is a very neat way to accommodate more than one graphic on one slide. A sample command is included for reference.

\newenvironment{wideitemize}{\itemize\addtolength{\itemsep}{10pt}}{\enditemize}

% Typically, the itemize enviroment does not provide adequate space between bullet points

% and hence can look a bit cluttered at times. The new "wideitemize" problem solves this issue.

% To use it just put your regular "\item" between "\begin{wideitemize}" and "\end{wideitemize}"

\hypersetup{% This can be left blank if information is not available

pdfauthor = {Author's name},

pdftitle = {Title for the pdf file},

pdfsubject = {},

pdfkeywords = {},

pdfcreator = {LaTeX with hyperref package},

pdfproducer = {dvips + ps2pdf}}

\mode<presentation> % This is a representative list of themes available is Beamer. Shreyes prefers the Warsaw theme, I prefer Pittsburg; it's much neater.

Also, CambridgeUS is quite popular among academics.

{

%\usetheme{Warsaw}

%\usetheme{Rochester}

%\usetheme{Madrid}

%\usetheme{Pittsburgh}

%\usetheme{Antibes}

%\usetheme{Montpellier}

%\usetheme{Berkeley}

%\usetheme{PaloAlto}

%\usetheme{Goettingen}

%\usetheme{Marburg}

%\usetheme{Hannover}

%\usetheme{Berlin}

%\usetheme{Ilmenau}

%\usetheme{Dresden}

%\usetheme{Darmstadt}

%\usetheme{Frankfurt}

%\usetheme{Singapore}

%\usetheme{Szeged}

%\usetheme{Copenhagen}

%\usetheme{Malmoe}

\usetheme{CambridgeUS}

\setbeamercovered{transparent = 28}

}

\usefonttheme{professionalfonts}

\usecolortheme[rgb={0.01,0.18,0.42}]{structure}

\useinnertheme{rounded}

\title{Sample Presentation}

\author{Author} \institute{My Institute}

\setbeamertemplate{footline}[frame number]

% This shows the slide number at the bottom of the page. May or may not be required depending on the theme you use.

\beamertemplatenavigationsymbolsempty

% This disables the navigation button in the lower right corner of the generated pdf file.

\begin{document}

\begin{frame} % Just so we are clear, a frame here means a slide.

\titlepage

\end{frame}

\AtBeginSection[] % This will display the list of contents before the beginning of each section

{

\begin{frame}

\frametitle{Outline}

\tableofcontents[currentsection]

\end{frame}

}

\begin{frame} % You can use this if you just want the list of contents to be displayed once at the beginning.

\tableofcontents

\end{frame}

\section{Introduction}

\begin{frame}

\frametitle{A simple sample frame}

\begin{wideitemize}

\item Point 1

\pause

\item Point 2

\pause

\item Notice the space between the points.

\pause

\item Comparatively, more generous than the normal ``itemize" enviroment

\pause

\item Also, the sequential appearance of points is handeled with \textbackslash pause command.

\end{wideitemize}

\end{frame}

\section{Part I}

\begin{frame}{Frames you can replicate} % You can include the frame title here as well.

\begin{itemize}

\item 1

\item 2

\item 3

\item 4

\item 5

\end{itemize}

\end{frame}

\begin{frame}{Another one you can replicate} % This one uses wideitemize

\begin{wideitemize}

\item 1

\item 2

\item 3

\item 4

\item 5

\end{wideitemize}

\end{frame}

\begin{frame}{Nested environments}

\begin{wideitemize}

\item This is wideitemize environment

\item You can nest environments to formulate a

hierarchical structure

\item See below

\begin{itemize}

\item This is itemize

\pause

\item Notice the space between bullets

\pause

\item The bullets can be changed to number using the enumerate environment

\pause

\item \LaTeX is cool

\pause
\item I am glad you chose to give it a

thought

\pause

\item It will change the way you work

\end{itemize}

\end{wideitemize}

\end{frame}

\section{Part II}

\begin{frame}{subfigure example}

\begin{figure}

\centering

\subfigure[Figure 1]{

\includegraphics[scale=0.4]{name_of_figure1_file.pdf}

\label{model1}} % Try putting labels for every figure, table, and equation. It will become easier to navigate if the document is huge.

\subfigure[Figure 2]{

\includegraphics[scale=0.4]{name_of_figure2_file.pdf}

\label{model2}}

\end{figure}

\end{frame}

\begin{frame}{Math in \LaTeX}

\begin{enumerate}

\item Use the ``align" or ``gather" environment for math

\item They are neater

\item You can also use the ``equation" and ``eqnarray" envrionments

\item Also, notice that this is the enumerate environment

\item Hence the numbers instead of bullets

\end{enumerate}

\end{frame}

\centering

\frame[plane]{

\vspace{2cm}

{\huge Thank you} % Huge makes the text size, well, Huge

\vspace{3cm}

\begin{flushright}

Author name and contact

\end{flushright}

}

\end{document}

Project Euler: problem 2

2011-09-16T17:15:00.000+05:30

Each new term in the Fibonacci sequence is generated by adding the previous two terms. By starting with 1 and 2, the first 10 terms will be:

1, 2, 3, 5, 8, 13, 21, 34, 55, 89, ...

By considering the terms in the Fibonacci sequence whose values do not exceed four million, find the sum of the even-valued terms.

# Inititae a vector x with two values 1 and 2, the starting points for the Fibonacci series

x <- c(1,2)

length(x)

# Take an object "i", with a starting value of 1.
# This object will be used to as an index for the vector "x".
# We continue to add# the (n - 1)^th term and the (n - 2)^th term
# to get the n^th term.

# This process continues as long as an element of vector x with
# index value "i" just crosses the 4,000,000 mark.

i <- 1

while (x[i] < 4000000){i <- i + 1

x.index <- length(x)

x[x.index + 1] <- x[x.index] + x[x.index - 1]}

# Sum the even values of the Fibonacci series thus obtained

sum(x[x %% 2 == 0])

Answer : 4613732

Project Euler: problem 1

2011-09-15T12:44:00.000+05:30

To be fairly honest (assuming there are degrees of honesty), I do know a little about math and programming but I don't know much math or any programming. I've loved math for a long time, but started to learn and understand fairly recently. So during the process of learning and understanding math and a little bit of programming, Shreyes and I thought of sharing our procedures. Most of these procedures aren't the most elegant solutions and at times are plain clumsy and inefficient. We plan to improve our programming skills as we go along.

We have recently started with Project Euler problems and will be posting some of the methods that we have used to arrive at a solution for each of the problems. We know only one language, R and hence our solutions are written in R.

So let's start with problem 1.

If we list all the natural numbers below 10 that are multiples of 3 or 5, we get 3, 5, 6 and 9. The sum of these multiples is 23.

Find the sum of all the multiples of 3 or 5 below 1000.

I tried two approaches to solve this problem and have detailed both of them below. By the time I was through with the problem, I found another approach.

First, let's see the one that I thought was quite efficient.

# Approach 1

# Create an empty vector of desired length where the results of the division will be stored

x.1 <- rep(NA, 999)

#NA is replicated 999 times

# For each integer from 1 to 999 (we want numbers below 1000), divide the integer # by either 3 or 5 and take the modulus.

# If the integer, say "i" is completely divisible by either 3 or 5, assign that value to the ith element of vector x; otherwise assign i'th value of x to equal to zero

for (i in 1:999){

if (i %% 3 == 0 || i %% 5 == 0) {x.1[i] <- i} else {x.1[i] <- 0}}

# "%%" is the modulus function, "||" is the symbol for the "or" command, x.1[i] calls the ith element of vector x.1

# Take all the non-zero values of x, i.e. those values for which the integer was perfectly divisible by either 3 or 5 and assign it to a separate vector

y.1 <- x.1[x.1 != 0]

# This assigns all the nonzero rows of vector x.1 to a vector y.1

# Take a summation of these values and this is the desired output.

sum(y.1)

# "sum" finds the sum of all the elements of vector y.1

Answer: 233168

Also, just to make sure that there is no confusion, the ".1" in x.1 and y.1 does not denote anything special, it is merely an assigning convention indicating that the variables were created for the first approach.

# Approach 2

# Take numbers from 1 to 333 and multiply each by 3 to get multiples of 3.

# We only take numbers till 333 because we have to find the sum of multiples of 3 (and 5) that are less than 1000.

x.2 <- 0

for (i in 1:333){x.2[i] <- i*3}

head(x.2) 

# head shows the first few rows of the object

# For multiples of 5, we take numbers till 199 and the last multiple of 5 below 1000 comes to be 995

y.2 <- 0

for (j in 1:199){y.2[j] <- j*5}

head(y.2)

# 15 is a common multiple of 3 and 5, and hence it will get included twice - once when we add the numbers that are multiples of 3 and once when we add numbers that are multiples of 5. So we need to subtract these 15 and its multiples.

z.2 <- 0

for (k in 1:66){z.2[k] <- k*15}

head(z.2)

# Final sum

sum(x.2) + sum(y.2) - sum(z.2)

Answer: 233168

There is another approach, which uses the funda of arithmetic progressions. Let's see if you can figure that out.

Project Euler: problem 5

2011-06-30T20:12:00.003+05:30

Calculate LCM of 'n' consecutive natural numbers using R
What is the smallest positive number that is evenly divisible by all of the numbers from 1 to 20?

Well I shall hit the nail right on the head and not beat around the bush. I am taking programming lessons on R from my pro bro(Utkarsh Upadhyay) who agreed on teaching me only if I would disseminate my learning(a paranoia all the open-source advocates share). Hence I shall populate the web with another link, which might help other dumb programmers like me.

Question:What is the smallest positive number that is evenly divisible by all of the numbers from 1 to 20?(Which is essentially the LCM of all the numbers from 1 to 20)

I came across this problem http://projecteuler.net/index.php?section=problems&id=5 here. You can choose any programming language to solve the problems given. I choose R.

My first attempt to solve this:

a <- 20

c <- 0
while ( c < 20){   // this loop stops when the value of ‘c’ becomes 20, i.e ‘a’ is divisible by all the numbers from 1 to 20
            c <- 0
            n <- 1
            a <- a + 1
            while ( n < 21 ){   // check if ‘a’ is divisible by all numbers from 1 to 20
                        if (a%%n == 0) c <- c + 1
                        n <- n+1
            }
}
print (a)

Brief note on the program above:

What I am essentially doing is checking if the number in the variable ‘a’ is divisible by all the numbers from 1 to 20, if its not then I am incrementing the value of ‘a’ and proceeding again with the loop(first while loop). So I would be checking for all the numbers starting 21 whether they are divisible by all the numbers from 1 to 20 and the program runs till ‘a’ takes the desired value.(which came out to be 232,792,560, after the computation was over).

This was a conservative way of getting the job done. It however took 9 hours of computation to blurt out the answer. Hmmm, well I could live with that number but just out of curiosity I asked Utkarsh if there was anything that I was missing. I just wish I were there to see the expression on his face, it would probably have been that of despair, or could also have been a hysterical laughter, I would never get to know that(sigh, Schrodinger's cat) but nevertheless lets focus on the task at hand. The suggestion Utkarsh gave was to use "recursive functions".

Revised program using recursive function and pro bro's help:

We essentially define 2 function and call one in the other. Its easier when you look at the code:

Defining a function lcm(a,b) and storing the codes in a file "LCM.R":

lcm <- function(a, b) {

if(a > b) { # Swap the numbers to keep the smallest number in ‘a’

a <- a + b

b <- a - b

a <- a - b

}

i <- 2

comb <- 1

while(i <= a) {

if(a %% i == 0 && b %% i == 0) {

# Accout for all the common factors

a <- a / i

b <- b / i

# Count common factors only once.

comb <- comb * i

}

else {

# i is not a common factor, carry on to the next number

i <- i + 1

}

return (comb * a * b) # For the non common factors, count all of them

}

A brief note on the above program:

What we have done here is we have defined a function lcm(a,b) as per our convenience and defined it such that it returns the LCM of ‘a’ and ‘b’. The logic used to calculate the LCM is what most of us have already used in class 5^th. Identify the common factors(which would be contained in ‘i’) and then to compute the LCM just use “comb <- comb * i”. Note that whenever I come across a common factor I am dividing both ‘a’ and ‘b’ by ‘i’ thus the values of ‘a’ and ‘b’ left in the end of the loop would be co-prime.(Think about this.!!). Therefor I am returning (comb * a * b), which would return the LCM of ‘a’ and ‘b’.

This program would be stored in a file let’s say “LCM.R”. Now whenever I have to refer to this function lcm(a,b) all I need to do is source this file “LCM.R” and I can conveniently use the function lcm(a,b) to get the LCM of ‘a’ and ‘b’. It would be as if the function lcm(a,b) always existed.

Now I can address a question for the novice programmers. Where do the 'a' and 'b' come from?

So whenever I use this self defined function lcm(a,b) I will use it in a program right.? so if I write

source('LCM.R') // this would allow you to use the function you defined in "LCM.R"

l <- lcm(6,8) // the value of 'a' would be 6 and 'b' would be 8

print (l)

I would get 24.

Defining another function lcm1(list.num) and storing the codes in "LCM2.R"

Now we come to the tricky part. What we will do now is use the function lcm(a,b), that we defined, and use it to compute LCM of a list of numbers.

source('LCM') //calling the file that stores the function lcm(a,b)

lcm1 <- function(list.num) {

LCM.so.far <- 1

for(next.number in list.num) {

LCM.so.far <- lcm(LCM.so.far, next.number) // Here lies the beauty

}

return (LCM.so.far)

}

A brief about the above program:

What we have done above is defined another function lcm1(list.num) which will take a list of numbers and blurt out the LCM.(Which is exactly what we want.!!). If we look at the ‘for’ loop defined we are running the loop for all the values in the list of numbers. Now the beautiful logic is in the line “LCM.so.far <- lcm(LCM.so.far, next.number)”, we have cleverly used the function defined earlier lcm(a,b) here. LCM.so.far would keep on updating it self as the loop runs with the next number in the list. Finally this function returns the LCM of the list of numbers that would be stored in ‘LCM.so.far’ at the end of the loop.(Think why.!!)

Now this function lcm1(list.num) and its definition would be stored in another file say “LCM2.R”. Similar to how we used the file “LCM.R” to use the function lcm(a,b) we can now source “LCM2.R” to use the function lcm1(list.num) that we defined.!!

So basically we have 2 function defined in 2 different files. To use the functions wll we need to do is source the files they are stored in.

Main Program:

source('LCM2') // this will call both the files(think why.!!)

l <- lcm1(1:20)

print (l)

Here we have sourced the file “LCM2.R” which would automatically open “LCM.R” too, since “LCM2.R” has to use lcm(a,b) defined in “LCM.R”(hope you catch the drift). Now we stored the LCM of the list of numbers from 1 to 20 (1:20) in ‘l’ and displayed ‘l’.

and TADA.!!!

The approximate computation time is <1 sec. And also we get a flexible functionality to compute the LCM of any consecutive list of natural numbers however long(not literally.!!)

Even if I consider that the computation took 1 sec, the program that I came up with took 9*60*60= 32,400 secs. Therefore the approximate efficiency enhancement achieved via this transition is 32,39,900% which is not bad I say..:-)

Critiques, abuses, banters, blessings are welcome.

P.S: Pardon me for the poor presentation and grammatical errors if any.