The Overlap Coefficient

I'm currently working with a client who is researching best covariance matrices for financial time series.  Specifically, they are looking at what best describes 1 month out of sample distributions.  They are not concerned with the means, just the variance.

A paper on the subject from the Journal of Portfolio Management -- "A Test of Covariance Matrix Forecasting Methods" by Zakamulin.

The test boils down to the sum of the squares of the differences in the upper half of the covariance matrices.  This is really just the Frobenius Norm of the difference between the matrices.

My issue with the Frobenius Norm has always been it is not, directly, a measure of the similarity between two distributions.  If we are using an estimated matrix to make investing decisions based on the distribution of likely outcomes, shouldn't we concern ourselves with the degree to which our estimated matrix can predict the distribution of outcomes?

Some thinking, noodling with code, thinking some more, and then a computer forced reboot (hosing said code) got me to actually Google the subject.  First thing that came up was this Cross Validated post.  Also this helpful presentation.  Seems what I was doing had been invented before (like that's a surprise).

"The Overlapping Coefficient (OVL) refers to the area under the two probability density functions simultaneously."


(Graph pulled from that Cross Validated post)

OVL is defined on the range of (0,1].  1 is a perfect fit with lower numbers a worse fit.

So how does the OVL compare to the Frobenius Norm (FN)?  When will their preferences in matrices differ?

We will make use of the R code from Cross Validated and use a univariate example as it is easy to visualize.  We will compare two distributions to a N(0,0.01) variable.  In the graph, the N(0,0.01) value is the solid line.  The distribution to compare is the dotted.

Full code at the end.

In the univariate case, FN is

FN is on the range of [0,∞).  0 is a perfect fit with higher numbers a worse fit.

Case 1:

X1 ~ N(0,0.01) vs X2 ~ N(0,0.0225)  -- Standard Deviation of 0.1 vs 0.15.

OVL = 0.8064

FN = 0.0125

Case 2:
X1 ~ N(0,0.01) vs X2 ~ N(0,0.0025)  -- Standard Deviation of 0.1 vs 0.05.

OVL = 0.6773

FN = 0.0075

The two competing distributions are both 0.05 different from the true standard deviation.  The OVL statistic prefers the Case 1 contender.  FN prefers the Case 2.  

Will the nature of the statistics cause the OVL to prefer higher SD values with the FN to prefer lower SD values?

Text in RED is the Standard Deviation being compared.

In general, we see the trend from upper left to lower right.  These statistics correlate on how they rank distributions in the way we expect.  The will both always pick the same preference on adjacent points.

We see that FN punishes over estimation of the SD.  That is, the higher numbers in red extend in a flatter arch down and to the right.

We see that OVL punishes under estimation of the SD.  That is the lower numbers in red extend in a steep arch down and to the right.

The FN is approximately indifferent between SD=0.02 and SD=0.14.  OVL would pick the 0.14 value.

The OVL is approximately indifferent between SD=0.07 and SD=0.14.  FN would pick the 0.07 value.

This leaves me unable to say which is better.  Intuitively, I prefer the OVL.  In a later post, I will develop a multivariate version of OVL using numerical integration.

As promised, R code for the above (yet again, hat tip to Cross Validated for a large portion of this):

min.f1f2 <- function(x, mu1, mu2, sd1, sd2) {

  f1 <- dnorm(x, mean=mu1, sd=sd1)

  f2 <- dnorm(x, mean=mu2, sd=sd2)

  pmin(f1, f2)

}

#CASE 1

mu1 <- 0

mu2 <- 0

sd1 <- .1

sd2 <- .15

xs <- seq(min(mu1 - 3*sd1, mu2 - 3*sd2), max(mu1 + 3*sd1, mu2 + 3*sd2), .01)

f1 <- dnorm(xs, mean=mu1, sd=sd1)

f2 <- dnorm(xs, mean=mu2, sd=sd2)

plot(xs, f1, type="l", ylim=c(0, max(f1,f2)), ylab="density")

lines(xs, f2, lty="dotted")

ys <- min.f1f2(xs, mu1=mu1, mu2=mu2, sd1=sd1, sd2=sd2)

xs <- c(xs, xs[1])

ys <- c(ys, ys[1])

polygon(xs, ys, col="gray")

print(paste("OVL:",integrate(min.f1f2, -Inf, Inf, mu1=mu1, mu2=mu2, sd1=sd1, sd2=sd2)$value,sep=" "))

print(paste("FN:",abs(sd1^2 - sd2^2),sep=" "))

#CASE 2

mu1 <- 0

mu2 <- 0

sd1 <- .1

sd2 <- .05

xs <- seq(min(mu1 - 3*sd1, mu2 - 3*sd2), max(mu1 + 3*sd1, mu2 + 3*sd2), .01)

f1 <- dnorm(xs, mean=mu1, sd=sd1)

f2 <- dnorm(xs, mean=mu2, sd=sd2)

plot(xs, f1, type="l", ylim=c(0, max(f1,f2)), ylab="density")

lines(xs, f2, lty="dotted")

ys <- min.f1f2(xs, mu1=mu1, mu2=mu2, sd1=sd1, sd2=sd2)

xs <- c(xs, xs[1])

ys <- c(ys, ys[1])

polygon(xs, ys, col="gray")

print(paste("OVL:",integrate(min.f1f2, -Inf, Inf, mu1=mu1, mu2=mu2, sd1=sd1, sd2=sd2)$value,sep=" "))

print(paste("FN:",abs(sd1^2 - sd2^2),sep=" "))

#Compare over range of SD values

x = matrix(seq(from=.01,to=.2,by=.01))

fn = apply(x,1,function(s1) abs(.01 - s1^2))

ovl = apply(x,1,function(s1) integrate(min.f1f2, -Inf, Inf, mu1=0, mu2=0, sd1=s1, sd2=0.1)$value)

toGraph = cbind(x,fn,ovl)

#indiffernce compared to SD=0.14

fn = abs(0.01 - .14^2)

ovl = integrate(min.f1f2, -Inf, Inf, mu1=0, mu2=0, sd1=.14, sd2=0.1)$value

plot(toGraph[,2],toGraph[,3],xlab="FN",ylab="OVL",ylim=c(0,1.1))

abline(h=ovl,lty="dotted")

abline(v=fn,lty="dotted")

text(toGraph[,2],toGraph[,3],toGraph[,1],col="red",pos=1)

Merging Data -- SAS, R, and Python

On analyticbridge, the question was posed about moving an inner join from Excel (which was taking many minutes via VLOOKUP()) to some other package.  The question asked what types of performance can be expected in other systems.  Of the list given, I have varying degrees of experience in SAS, R, and Python.

In the question, the user has about 80,000 records that match between 2 tables and need to be merged by a character key.

To start, let's create 2 tables in SAS.  Each will have 150,000 records and there will be 80,000 overlapping records.  We will randomize the tables (obviously the merge is trivial when the tables are sorted with the matching records on top) to be more "real world."

data A;

format address $32. a $8.;

a = "1";

do i=1 to 150000;

       rand = ranuni(123);

       if i <=80000 then do;

              address = put(i,z32.);

       end;

       else do;

              address = "AAA" || put(i,z29.);

       end;

       output;

end;

run;

data B;

format address $32. b $8.;

b = "1";

do i=1 to 150000;

       rand = ranuni(234);

       if i <=80000 then do;

              address = put(i,z32.);

       end;

       else do;

              address = "BBB" || put(i,z29.);

       end;

       output;

end;

run;

proc sort data=a;

by rand;

run;

proc sort data=b;

by rand;

run;

Now SAS pages tables to the hard drive.  Good and bad.  Python and R will be starting with the tables in memory, so we use SASFILE to load the tables into the main memory.  Also note that SAS is writing the result table to the HD.  We will do the same in Python and R.

594  sasfile a load;

NOTE: The file WORK.A.DATA has been loaded into memory by the SASFILE statement.

595  sasfile b load;

NOTE: The file WORK.B.DATA has been loaded into memory by the SASFILE statement.

596

597  proc sql noprint;

598  create table temp.tableA_B as

599  select a.address,

600         a.a,

601         b.b

602      from a inner join

603           b

604        on a.address = b.address;

NOTE: Table TEMP.TABLEA_B created, with 80000 rows and 3 columns.

605  quit;

NOTE: PROCEDURE SQL used (Total process time):

      real time           0.09 seconds

      cpu time            0.09 seconds

So SAS took 0.09 seconds.  Much faster than the many minutes in Excel.

In R this is trivial.  I wrote the tables to csv files (not shown).  So we will read them in, do the merge, and then save the result as an .rda file.

> system.time({

+   a = read.delim('c://temp//a.csv',header=T,sep=',')

+   b = read.delim('c://temp//b.csv',header=T,sep=',')

+ })

   user  system elapsed

   9.03    0.06    9.09

> system.time({

+   m = merge(a[c("address","a")],b[c("address","b")])

+ })

   user  system elapsed

   2.15    0.02    2.17

> system.time({

+   save(m, file="c://temp//tableA_B.rda")

+ })

   user  system elapsed

   0.21    0.00    0.20 

R took 2.17 seconds for the merge and 0.20 seconds to write.  A total of 2.37 seconds.

I am least familiar with the optimal way to do this in Python.  I have a question to my Python guru about the optimal way to do the merge.  For the time being, my attempt is here.  The basics is to read the files into a Dictionary with the address string as the key and a basic object as the value.  Then iterate over the keys in 1 table and see if they are in the second table.  If so, add the merged data to a new dictionary.

Outputting the merged data might be faster in a list instead of a dictionary.  The index hash is not being built in SAS or R.

C:\Temp>c:\Python27\python.exe merge.py

2.4200000763 Starting Merge

Took 0.2660000324 seconds to merge, 1.1210000515 total with pickle

Python took 0.27 seconds for the merge and a total of 1.12 seconds for the merge and write.

A final note, here is a screen grab of the resulting file sizes.  R wins hands down -- I assume there is some compression going on there.

This example is fairly trivial.  Hopefully, someone will find it useful while trying to learn one of these languages.

Expected Shortfall Portfolio Optimization in R using nloptr

I have previously done examples of QP optimization in for financial portfolios.  I am not a huge fan of variance optimization in finance.  Return distributions are not normal, are often skewed, and are usually leptokurtic.  In plain speak, the distributions have fat tails and while the mean may be 0, the median is shifted to one side or the other.

Because of the shape of return distributions, variance optimization fails to capture the true risk of a portfolio. The risk measure I prefer to use is Expected Shortfall.  Expected Shortfall is defined as the mean of the left tail of the distribution, below some given alpha.

Beyond VaR:

We call the value at the apha Value-at-Risk (VaR).  VaR has been much maligned over the last few years -- not always for the right reasons.  Some people confuse all VaR measures as assuming a normal distribution. Not so.  The definition is that is it some left tail alpha value of the expected P/L distribution.  No where does it assume a normal distribution -- people just wrongly use a normal distribution when calculating VaR.

Jorion's definition of VaR fromValue at Risk: The New Benchmark for Managing Financial Risk, 3rd Edition says (I'm paraphrasing) VaR is the maximum expected loss at a given level of certainty during normal market conditions.  That often is misinterpreted by management as THE number of the maximum loss.  Jorion's definition has a positive spin that pushes the misinterpretation.

I prefer to but a negative spin on it.  VaR is the minimum expected loss when you have a bad day.  If I have a day that is, at least, in the alpha percential of crappy-ness, then I will lose MORE than VaR.

VaR is the floor, not the ceiling.  That's where it has been misused.

So why not optimize on VaR?  Short answer is that VaR is highly nonlinear.  Add nonlinear instruments to your portfolio and you can no longer guarantee the concavity of the objective function.

Enter ES:
Expected Shortfall (ES) is the average of the values above the VaR value.  That is, when I'm having a really bad day, what is my expected loss.  If I say a 1 in 20 day is bad (alpha=.05) then I should have 12-13 bad days a year.  If my portfolio was static and the return distributions unchanging (not likely), then the average of those 12-13 days should center on ES.

ES is still not perfect because you don't know how far to the left you can go.  But it does factor in the fat, skewed tail that we see.  It is a good starting point for risk analysis.

Because ES is a mean, it is guaranteed to be a concave function (I'm still looking for the reference -- I've seen it before and will update the post if I find it).  So ES is a better measure of risk than VaR and we can trust optimization results from it.  So let's build an optimization routine for ES.

nloptr:
nloptr is a nonlinear optimization library in R wrapping the GNU NLopt library.  ES, while well behaved, is nonlinear.  nloptr provides a number of nonlinear solvers.  It's what we want.

For our test case, we will simulate a 4 variable normal distribution with 10,000 draws (correlation given below).  We will build the ES function and a gradient function for ES.  We will constrain portfolio weights to be between [0,1].  We will impose an equality constraint that the sum of the weights = 1.

require(MASS)

require(nloptr)

#Covariance structure for the simulation

#give everything a std=.1

n = 4

corr = c(1, .7,  .5,  .1,

              .7,  1,  .4, -.1,

              .5, .4,   1,  .6,

              .1, -.1, .6,  1)

corr = matrix(corr,n,n)

std = matrix(0,n,n)

for(i in 1:n){      

       std[i,i] = .1

}

cov = std %*% corr %*% std

#Simulate 10,000 draws

sim = mvrnorm(n=10000,rep(0,n),cov)

#feasible starting values of equal weights

w = rep(1/n,n)

#ES function.  Mean of values above alpha

es = function(w,sim=NA,alpha=.05){

       ret = sort(sim %*% w)

      

       n = length(ret)

       i = alpha * n

      

       es = mean(ret[1:i])

      

       return(-es)  

}

#linear equality constraint

#note: nloptr requires all functions to have the same signature

eval_g0 <- function(w,sim=NA,alpha=NA) {

       return( sum(w) - 1 )

}

#numerical approximation of the gradient

des = function(w,sim=NA,alpha=.05){

       n = length(w)

       out = w;

       for (i in 0:n){

              up = w;

              dn = w;

              up[i] = up[i]+.0001

              dn[i] = dn[i]-.0001

              out[i] = (es(up,sim=sim,alpha=alpha) - es(dn,sim=sim,alpha=alpha))/.0002

       }

       return(out)

}

#use nloptr to check out gradient

check.derivatives(w,es,des,sim=sim, alpha=.05)

#function to optimize -- a list of objective and gradient

toOpt = function(w,sim=NA,alpha=.05){

       list(objective=es(w,sim=sim,alpha=alpha),gradient=des(w,sim=sim,alpha=alpha))    

}

#equality constraint function.  The jacobian is 1 for all variables

eqCon = function(w,sim=NA,alpha=.05){

       list(constraints=eval_g0(w,sim=NA,alpha=.05),jacobian=rep(1,4))     

}

#optimization options

opts <- list( "algorithm" = "NLOPT_LD_SLSQP",

              "xtol_rel" = 1.0e-7,

              "maxeval" = 1000)

#run optimization and print results

nl = nloptr(w,toOpt,

              lb = rep(0,4),

              ub = rep(1,4),

              eval_g_eq=eqCon,

              opts=opts,

              sim=sim,alpha=.05)

print(nl)

s = nl$solution

obj = nl$objective

To confirm the results, I ran this 100 times and averaged the resulting weights and objective value.  I did the same in SAS IML.  UPDATE: SAS code located here.

From R we get

> apply(s,2,mean)

[1] 0.101697131193604 0.420428046895474 0.000000000004212 0.477874821906712

> mean(obj)

[1] 0.1371

In SAS we get

wgt

0.1015076 0.4218076 -4.75E-20 0.4766848

obj

0.1375891

I am comfortable with the results.

Final Thoughts:
I was only able to get the SLSQP routine to converge.  This routine uses a local quadratic approximation of the function, runs a QP to update the variables, and then iterates.  I am not sure why the others failed.

The number of routines available for problems with equality constraints are limited.  On top of that, the NLopt documentation gives a number of other routines that should accept linear constraints but the R implementation throws an error.  A number of augmented LM routines are available for equality constraints, but again the only one that worked used SLSQP as the sub-optimizer and produced the same result as SLSQP (in about 5x the time).

The time to run the optimization in R is high.  During the 100 iteration sample it took from 1.5-10 seconds per loop depending on the internal iterations needed by nloptr.  In all it took about 10 minutes.  The SAS code ran in 45 seconds.  I see the following as reasons:

1. SLSQP is not as efficient as other routines.  It has to approximate a hessian matrix numerically.  Finding another routine that reliably converges could be a big help.
2. My gradient function is numeric.  Each time the gradient is computed it takes a large number of FLOPs.  This is the same in SAS where the IML optimizer calculates numeric derivatives.
3. As previously discussed, the linear algebra matrix multiplication in R is slow.  This code relies heavily on it.  SAS IML is optimized for matrix math.

If anyone has ideas on how to speed up the R processing, I would love to hear them.  Maybe there are other, more efficient, nonlinear optimization libraries in R.  I'm sure there are parts of my code that can be sped up as well. 

Performance with foreach, doSNOW, and snowfall

Is it just me, or does the performance of the foreach package with a doSNOW backend operating on a socket grid suck?

Here at work, I am helping to setup a cluster of Windows machines for distributed R processing.  We have lots of researchers running code that takes hours to complete and are essentially large for loops with lots of analysis in between.  These guys and gals are not hard core programmers, so there is lots of interest in foreach (as opposed to something like RMPI).

I have successfully setup a POC grid between mutliple machines using sockets and public key authentication.  Assuming we use this, I'll post a how-to, as there is not much on the web on how to get it working on Windows.

In the meantime, I am testing performance.  There is something going on with foreach that I do not understand.  Performance numbers are really bad.

Can anyone explain what is going on here?

> require(doSNOW)
Loading required package: doSNOW
Loading required package: foreach
foreach: simple, scalable parallel programming from Revolution Analytics
Use Revolution R for scalability, fault tolerance and more.
http://www.revolutionanalytics.com
Loading required package: iterators
Loading required package: snow
> require(snowfall)
Loading required package: snowfall
>
> sfInit(parallel=TRUE,socketHosts=rep("localhost",3))
R Version:  R version 2.15.0 (2012-03-30)
snowfall 1.84 initialized (using snow 0.3-9): parallel execution on 3 CPUs.
> cl = sfGetCluster()
>
> f = function(x) {
+    sum = 0
+    for (i in seq(1,x)) sum = sum + i
+    return(sum)
+ }
>
> registerDoSNOW(cl)
>
> out = vector("logical",length=10000)
> system.time( (for (i in seq(1,10000)) out[i]=f(i) ))
   user  system elapsed
  25.99    0.00   25.99
>
> system.time( (out = lapply(seq(1,10000),f) ))
   user  system elapsed
  26.55    0.00   26.55
>
> system.time( (out = parLapply(cl,seq(1,10000),f) ))
   user  system elapsed
   0.02    0.00   15.85
>
> system.time( (out = foreach(i=seq(1,10000)) %dopar% f(i) ))
   user  system elapsed
   6.64    0.42   98.31
>
> getDoParWorkers()
[1] 3

EDIT: HA!  Figured it out.  foreach is not very efficient in communicating tasks as compared to par*apply().  The time to communicate the process overwhelmed the actual processing time.

When I change the code to this, it runs fast (about the same as parLapply()):

> system.time( (out = foreach(i=seq(0,9),.combine='c') %dopar% {
+    apply(as.array(seq(i*1000+1,(i+1)*1000)),1,f)
+ }))
   user  system elapsed
   0.00    0.00   14.03

Calling R from SAS IML Studio

I am playing around with SAS IML Studio 3.4.  For those that do not know, IML (Interactive Matrix Language) is the Matlab-esk language from SAS.  It opperates from normal SAS code through the PROC IML procedure.  A new (to me at least) UI has been developed for analysts called IML Studio.  IML Studio uses a superset of the IML language called IMLPlus.  I'll be digging into it (and the goodies like linked graphs, Java integration, and the ability to call 3rd party dll's) later.

One of the more recent additions to the IML and IMLPlus languages is the ability to run SAS routines from within IML.  At the same time this functionality was added, SAS also added the ability to call R from within IML.  You can now pass IML matrices back and forth betwen R matrices and SAS Datasets back and forth to R Data Frames (and other types).

Having never done this, I fired up IML Studio and set out to learning.

First, save the macros created in my last post into the an Autocall library.  You can modify the autocall libraries by modifying the sasv9.cfg file and adding the path to the SASAUTOS list.  Mine looks like this:

-SET SASAUTOS (
       
        "C:\Users\pazzula\Documents\My SAS Files(32)\9.3\macros"
        "!SASROOT\core\sasmacro"
        "!SASROOT\accelmva\sasmacro"
        "!SASROOT\dmscore\sasmacro"
        "!SASROOT\ets\sasmacro"
        "!SASROOT\iml\sasmacro"
        "!SASROOT\stat\sasmacro"
        )

The SAS file for the macro can be downloaded here.

To submit SAS code, surround the code with "submit;" and "endsubmit;".  This piece will download data for the SPY ETF:

submit;

%get_stocks(spy,25MAY2010,);

endsubmit;

Next, let's create 2 vectors, X and Y.  Make Y and linear function of X.

x = (1:10)`;

y = 1 + 3*x;

e = j(10,1,0);

do i=1 to 10;

               e[i] = .5*rannor(12345);

end;

y = y + e;

 Nothing hard about that. Those new to IML will want to know that ` is the transpose operator and "j(n,m,value)" creates a matrix (n x m) filled with "value."

Exporting IML matrices and SAS Data Sets to R is straight forward.  Use the modules ExportMatrixToR() and ExportDataSetToR().

run ExportMatrixToR(x,"x");

run ExportMatrixToR(y,"y");

run ExportDataSetToR("returns","returns");

The second parameter to each module is the name to give the object in R.  To call R, we again use "submit" and "endsubmit," only this time we add "/ R" to the submit line.  So let's run a linear model on y~x, create an XTS object from the returns Data Frame, chart the cumulative returns of SPY and create an AnnualizedReturn table.

submit /R;

require(xts);

require(PerformanceAnalytics);

m = lm(y~x);

summary(m);

returns = xts(returns$spy,returns$Date);

colnames(returns) = {"SPY"};

chart.CumReturns(returns[,"SPY"],main="Total Return");

table.AnnualizedReturns(returns);

endsubmit;

Produces

Call:
lm(formula = y ~ x)
Residuals:
Min 1Q Median 3Q Max
-0.79782 -0.04944 0.04503 0.17198 0.33329
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.85948 0.22469 3.825 0.00505 **
x 3.02539 0.03621 83.548 4.7e-13 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.3289 on 8 degrees of freedom
Multiple R-squared: 0.9989, Adjusted R-squared: 0.9987
F-statistic: 6980 on 1 and 8 DF, p-value: 4.698e-13
                                                SPY
Annualized Return                0.1041
Annualized Std Dev              0.1956
Annualized Sharpe (Rf=0%) 0.5323
and:

The SAS Data Set contained a column called Date that had a SAS Date format applied.  During the conversion to R Data Frame, SAS was nice enough to convert that column into an R date.  

That's pretty much it.  It's pretty straight forward.  Personally, I'm excited about this.  There are some things, like data manipulation, that SAS is way better than R at.  But then there are things that R gives me that I have to work to code in SAS (like easy functions for portfolio analytics).  Now I get the best of both worlds.