Showing posts with label R. Show all posts
Showing posts with label R. Show all posts

Wednesday, July 24, 2013

Casting and Melting with Paired Data

Today's post is not about economics, rather it's a note from an R programming struggle that may be helpful for fellow undergraduate researchers.

I'm often testing forecasting models, and what this ends up creating is a bunch of "forecasted" variables that are paired with the "actual" values. R has fabulous faceting capabilities, and I have often wanted to reshape the data in a way where the category of forecasted variable as an identifier, and then two columns that list the forecasted and actual variables. In other words, if the code starts from something like


       aAct      aPred       bAct      bPred id
1 1.2076384 -0.6735547  1.4994464 -1.0691975  1
2 0.4999706 -0.7188215 -0.3601551  0.7224729  2
3 1.0340859 -0.1108304 -0.5941295  0.5027085  3

And I want to convert it where one column has an id, another one identifies whether I'm forecasting a or b, and a third column that has the forecasted value, and then a fourth column with the actual value.

The procedure in R involves "melting" the data frame and then "casting" it. Melting is rather simple -- you provide a set of identifiers, and then the data frame is melted down to only that identifier, the values, and another indicator variable that tells you what the value is supposed to represent. In the above example, if we let df be the data frame described above, I would run:


df.m = melt(df, id.vars = 'id')

   id variable      value
1   1     aAct  1.2076384
2   2     aAct  0.4999706
3   3     aAct  1.0340859
4   1    aPred -0.6735547
5   2    aPred -0.7188215
6   3    aPred -0.1108304
7   1     bAct  1.4994464
8   2     bAct -0.3601551
9   3     bAct -0.5941295
10  1    bPred -1.0691975
11  2    bPred  0.7224729
12  3    bPred  0.5027085


Now I need to "unmelt" part of the data frame to get the forecast/actual pairings. In R, this is known as casting and I know that I personally had a pretty hard time decoding the documentation. The function goes along as

cast(df.m, <IDENTIFIERS> ~ <VALUES>)

The second part is known as the casting formula and is the part that I have struggled with. But in its most simplest form, the casted frame will look like something with all the identifiers added together as uniquely identifying units ,and then the <VALUES> variables being the labels for the actual value column. If that sounded confusing, I apologize. Perhaps solving the example would help.

First, I need to find a way to identify whether a row is looking at a or b, and whether it is a forecast or an actual variable. So I first create these variables:

df.m$var = substring(df.m$variable, 1, 1)
df.m$type = substring(df.m$variable, 2) 

Which gives me the data frame:


> df.m
   id variable      value type var
1   1     aAct  1.2076384  Act   a
2   2     aAct  0.4999706  Act   a
3   3     aAct  1.0340859  Act   a
4   1    aPred -0.6735547 Pred   a
5   2    aPred -0.7188215 Pred   a
6   3    aPred -0.1108304 Pred   a
7   1     bAct  1.4994464  Act   b
8   2     bAct -0.3601551  Act   b
9   3     bAct -0.5941295  Act   b
10  1    bPred -1.0691975 Pred   b
11  2    bPred  0.7224729 Pred   b
12  3    bPred  0.5027085 Pred   b

Now I can cast the frame. In this case, I would use the formula

df.mc = cast(df.m, id + var ~ type)

This is how you interpret the formula. Id + var means that every observation is uniquely identified by it's id code and the variable we're forecasting -- a or b. Then "type" on the right side represents the new variable names that will be filled by the values.

Hope this is useful to others so they don't end up spending hours agonizing over the issue as did I.


Thursday, August 2, 2012

Inaccurate Estimation in a Gaussian World

How accurate is estimation in a Gaussian world?

Critics of finance lavish lots of attention on "fat-tail distributions" and how they call into question the way finance deals with low probability events. Nicholas Nassim Taleb is particularly angered by the way financiers use the Gaussian bell curve, affectionately known as the Great Intellectual Fraud, to "predict" and optimally hedge bets. While sympathetic to this argument, I still find the Gaussian bell curve to be a useful tool to help demonstrate how fragile and unpredictable low probability events are even in a normal world. Working on the problem also proves to be a convenient time to orient myself in R, a statistical package that I will likely be using in my undergraduate research at the University of Michigan this fall.

So here's the problem I want to look at:

Given a sample from a normally distributed population
1) How accurate of an estimate of the population standard deviation is the sample standard deviation?
2) As a result of (1), by how much do you over or under-estimate the probability of tail events?
3) How does the over or underestimation change as the event you're trying to estimate becomes more extreme (ie higher sigma event?)

I've previously written about (2), but the point of this post is to run the full simulation and what are the results.

To start, I generate a 500 element population, which has both a standard deviation and mean of 1.


From here, I take 1000 samples of 50 each, and from each sample I calculate the standard deviation. Below is the distribution of the percent errors of those standard deviations estimates.

As you can see, in spite of the fact that the standard deviation was measured 1000 times, there's still a substantial amount of spread in the distribution of standard deviations. While they average a 3% underestimate, they range from an underestimation of around 40% to an overestimation of over 20%.

This spread in the standard deviation estimate is particularly worrisome when one starts estimating the probability of low probability tail events. Below is the ratio of the actual left tail probability of a 3 sigma event divided by the estimated probability based on the standard deviation estimates. A big number implies an underestimation of the tail risk, while a small number implies an overestimation of the tail risk. Note the skew.

So when we're talking about a 3 sigma event, there's a very sizable risk of a dramatic underestimate. While most of the underestimation is concentrated around 1.14, which corresponds to a 13% underestimation, the ratio can go up to around 15, suggesting that risks can be massively underestimated even in a Gaussian world. Astute readers may note that the estimation error distribution looks to be log-normal, which at least has finite variance. But if a normal population can lead to a log-normal distribution of estimation error, it means that they end up with even more potential for underestimation. As pictured below, a log-normal population means your average underestimate is even higher at around 5 times, while the tail is even fatter.

Inline image 1
Now that we've explored part 2, we can look at part 3; how does this change at different levels of standard deviation? I generate a population with a standard deviation , and collect data on it as per the sampling procedure above. I do this over and over again at different populations at different standard deviation levels. Surprisingly, I don't get a robust result for the estimation of error, as it's highly dependent on the population generated. However, most of the time, I do get the general trend for the mean of the log of the error ratio, pictured below. Note that it becomes more and more negative, suggesting that the mean underestimation goes down as the severeity of the tail event increases. This is likely because of the way the skewness works out with the over and under estimation of the standard deviation.
But is this necessarily good news? Not really. While the mean might be tending towards less underestimation, the proportion of underestimated risk stays relatively constant.

While the mean value of the underestimation given the risk is underestimated shoots up.

From these graphs, we can see that the average magnitude of the underestimation is increasing very quickly, while the proportion of underestimation stays relatively constant. It is actually increasing faster than it looks because the y-axis is the log of the ratio. Combined with the fact that the proportion of underestimation is staying constant, this implies the distribution is getting flatter and more dispersed, opening the possibility of catastrophic loss. The mean underestimation likely understates the actual damage that would result, as only one extremely bad underestimation can send the firm bankrupt, along with possibly the rest of the industry.

These graphs show with great detail at how insufficient risk measurement techniques like VaR or even ES are. Small probabilities are hard to estimate, even in a normal, Black-Swan free world. But add some grey swans, fragile balance sheets, and large banks, it's a time bomb that's just waiting to explode.