There are two popular parametric approaches to EVT that we will cover in this post. The first is called Block Maxima method. The second is called Peak-Over-Threshold (POT).

Block Maxima

Main idea:

The main idea of Block Maxima is that given we have some historical raw data, we break it into k blocks (e.g. calendar month/quarter/year), and look at the maximum value in each block and gather them as sample data to model tail risk distributions.

• The Gumbel

• The Fréchet

• The Weibull

The nice thing is that all three can be combined into a single family of continuous distributions known as Generalized Extreme Value (GEV) distribution:

Lets now have a look at what the three distributions look like.  In R package fExtremes has function dgev to plot the GEV density.  

require(fExtremes)
par(mfrow=c(1,3))
x <- seq(from = -7.5, to = 7.5, length.out = 550) 

# plot Weibull
plot(x, dgev(x, xi = 0.5, mu = 0, beta =1),
     type = 'l', lty = 1, lwd = 2,
     col = 'steelblue', ylab = 'Density')

# plot Gumbel
plot(x, dgev(x, xi = 0, mu = 0, beta =1),
     type = 'l', lty = 1, lwd = 2,
     col = 'steelblue', ylab = 'Density')

# plot Frechet
plot(x, dgev(x, xi = -0.5, mu = 0, beta =1),
     type = 'l', lty = 1, lwd = 2,
     col = 'steelblue', ylab = 'Density')

Below are the density plots:

At this point I hope the main idea of Block Maxima approach to EVT are well understood and I will move on to playing with data and seeing how we can put the theory to work.

Data Fitting – Block Maxima:

We now know that we want to fit a GEV distribution to some data set.  The issue is how to do it.  We are trying to model maxima so lets try the following.  Lets take some historical data, break it into k blocks, and look at the maximum value in each block and consider that to be our sampled value that we are trying to model.  Lets assume we are working with a financial time series and we are looking at daily returns.  Furthermore assume we are trying to model loss distribution so we will look at only negative returns on a daily basis.  Below is a stylized chart assuming looking at only losses (represented as positive values).  We are looking at 3 months worth of data.  We can break the data set into monthly blocks and take measure the maximum in each month.  We highlight the three observations below in red.  The three observations would be used to fit the GEV distribution. 

In our example we will be working with daily returns on SPY.  We can load n years worth of data and break it into monthly blocks.  So we will have k=n/12 blocks.  For this example we will use R’s quantmod package to load data for 13 years.  We then calculate maximum daily loss for each month and also keep the day of the observation even though it is not necessary.

require(quantmod)
require(xts)

SPY_prices <- Ad(getSymbols('SPY',
                            from = '2006-01-01',
                            to = '2019-01-01',
                            auto.assign = FALSE))

returns_xts <- diff(log(SPY_prices), k = 1)[-1]*100  # Daily log returns

# below reverses sign of returns_xts so loss is a positive number
# then collect max value for each month

monthly_max_daily_loss <- do.call(rbind,
                                  lapply(split(x = -1*returns_xts, 'months'),
                                        function(x) x[which.max(x)]))

Since we are trying to maximize the likelihood, but our optimization routines are for minimization routines, we will implement the negative log likelihood function and minimize it. Mathematically the two are equivalent.

# likelihood function
gev_loglik <- function(params, data){
    mu <- params[1]; sigma <- params[2]; xi <- params[3]
    data <- coredata(data)
    m <- nrow(data)
    
    if((sigma <=0) | (xi <= -1)) return(1e6)  # if outside bounds return large value
    term1 <- -m*log(sigma)
    term2 <- (data-mu)/sigma
    
    # case xi = 0 (ie Gumbel distribution)
    if(abs(xi)<= .0001){
        result <- term1 - sum(term2) - sum(exp(-term2))
    } else {
        if(any(1+xi*term2 <= 0)) return(1e6)
        result <- term1 - (1+(1/xi))*sum(log(1+xi*term2))-sum((1+xi*term2)^(-1/xi))
    }
    return(-result)
}

Now let’s try to minimize the function:

my_gev_fit <- optim(par = c(1,1,.05), fn = gev_loglik, method = "Nelder-Mead", data = monthly_max_daily_loss)
my_gev_fit

$par
1.26096477678288   0.799888376043898   0.275120760011372
$value
235.677220980881
$counts
function   70 
gradient   <NA>
$convergence
0
$message
NULL

We can compare this result with fExtremes’ gevFit function.

gevFit(coredata(monthly_max_daily_loss), block = 1, type = 'mle')

Title:
  GEV Parameter Estimation 
Call:
  gevFit(x = coredata(monthly_max_daily_loss), block = 1, type = "mle")
Estimation Type:
   gev mle 
Estimated Parameters:
        xi        mu      beta 
 0.2751779 1.2611064 0.7999340 

You can see that our estimated parameters are close.

Now that we have our estimated parameters, we can calculate various risk measures. For example, we can compute Value at Risk (VaR) and Expected Shortfall (ES).

The inverse of the cumulative distribution function is given as:

The implementation is below:

GEV_VAR <- function(params, alpha = .05){
    # to vectorize the function we need vectors of estimated coefficients
    mu    <- rep(params[1], length(alpha))
    sigma <- rep(params[2], length(alpha))
    xi    <- rep(params[3], length(alpha))
    
    y <- -log(1-alpha)
    result <- ifelse(abs(xi) < 0.0001, mu -sigma*log(-y), mu -sigma/xi*(1-y^-xi))
    return(result)
}

Now we can extract our fitted coefficients and feed them into our VaR function and estimate the VaR at various alphas. Please note that in this implementation, alpha of 5% means it is the 95th percentile of the distribution.

my_fit_vals <- my_gev_fit$par
GEV_VAR(my_fit_vals, alpha = c(.05, .025, .01))

4.93706812468633   6.34878189965057   8.66305966965268

Let’s check our results versus fExtremes’ native quantile function using the fitted values from the gevFit function:

# lets check the results verus fExtremes' quantile function
qgev(p = 1- c(.05, .025, .01), xi = 0.2751779, mu = 1.2611064, beta = 0.7999340)

4.93682630963963   6.3484731241556   8.66265699310054

Our results are very close.

Finally, let’s write a function to get the Expected Shortfall (ES). And let’s calculate the ES at 1% level. We will implement a numerical integration here.

GEV_ES <- function(params, alpha = .05){
    my_fun <- function(x){GEV_VAR(x, params = params)}
    result <- 1/alpha * integrate(my_fun, lower = .00001, upper = alpha, stop.on.error = FALSE)$value
    return(result)
}

GEV_ES(my_fit_vals, alpha = .01)

12.49453391079

Notice that at 1% the VaR is 8.66% but the expected shortfall (conditional VaR) is much more significant, 12.50%.

Peak-Over-Threshold

Main idea:

A more common approach to EVT in finance is the so-called Peak-Over-Threshold (POT) method. It must be obvious to the reader that we have been very wasteful with our data up to this point. We discarded all observations in every month except for one extreme.

Occasionally, we may observe multiple “extreme” observations in a given time block, and other months may be quiet. This may create instances where we under-represent the tail of the distribution when estimating the GEV model. For example, using Block Maxima method, we would’ve picked the losses circled in red (i.e. largest loss in each calendar month) to model tail distributions:

As shown in above graph, we would’ve ignored many extreme losses occurred in March, while using insignificant losses from January and February when we use Block Maxima method.

In POT we instead retain observations that exceed some threshold as shown below.

Data Fitting – POT:

A theorem in EVT called Pickands-Balkema-de Haan states that if Block Maxima have approximately GEV distribution, the values in excess of some threshold (for a large enough threshold) has a limiting distribution given as below:

This distribution is referred as the Generalized Pareto Distribution (GPD).

Like what we did for GEV, let’s plot the density functions using fExtremes package.

par(mfrow=c(1,3))
x <- seq(from = -7.5, to = 7.5, length.out = 550) 


plot(x, dgpd(x, xi = 0.5, mu = 0, beta =1),
     type = 'l', lty = 1, lwd = 2,
     col = 'steelblue', ylab = 'Density')

plot(x, dgpd(x, xi = 0, mu = 0, beta =1),
     type = 'l', lty = 1, lwd = 2,
     col = 'steelblue', ylab = 'Density')

plot(x, dgpd(x, xi = -0.5, mu = 0, beta =1),
     type = 'l', lty = 1, lwd = 2,
     col = 'steelblue', ylab = 'Density')

gdp_loglik <- function(params, threshold, data){
    xi <- params[1]; tau <- params[2]
    data <- data - threshold # y = x-u
    data <- coredata(data[data >0]) # only need data that exceeds threshold
    m <- nrow(data)
    
    if((tau <=0) | (xi <= -1)) return(1e6)  # if outside bounds return large value
    term1 <- -m*log(tau)
    
    #  case xi = 0
    if(abs(xi) < 0.0001){
        result <- term1 - 1/tau *sum(data)
    } else {
        if(any(1+xi*data/tau <=0)) return(1e6)
        result <- term1 - (1+1/xi)*sum(log(1+xi*data/tau))
    }
    return(-result)
}

Now, let’s load some SPY data again and fit the parameters.

require(quantmod)
require(xts)

SPY_prices <- Ad(getSymbols('SPY',
                            from = '2006-01-01',
                            to = '2019-01-01',
                            auto.assign = FALSE))

returns_xts <- -1*diff(log(SPY_prices), k = 1)[-1]*100  # Daily log returns
names(returns_xts)<- 'SPY_Returns'

Using a threshold of 0.85 we can run our optimization.

my_threshold <- .85
#my_gpd_fit
my_gpd_fit <- optim(par = c(.5,.5), 
                    fn = gdp_loglik, 
                    method = "Nelder-Mead", 
                    threshold = my_threshold, 
                    data = returns_xts)
my_gpd_fit

$par
0.150591120379977   0.861181165068411
$value
476.579589960172
$counts
function 55 
gradient <NA>
$convergence
0
$message
NULL

And let’s compare it to fExtremes package output.

gpdFit(coredata(returns_xts), u = my_threshold, type = 'mle')

Title:
   GPD Parameter Estimation 
Call:
  gpdFit(x = coredata(returns_xts), u = my_threshold, type = "mle")
Estimation Method:
   gpd mle 
Estimated Parameters:
        xi      beta 
 0.1507411 0.8610883

Again, we are close.

Now that we have our estimated parameters, we can move on to calculating our two risk measures: VaR and ES.

The quantile function is given as:

GPD_VAR <- function(params, data, threshold, alpha = .05){
    n <- nrow(data) 
    k <- sum(data > threshold)  # number of obs that exceed threshold
    # to vectorize the function we need vectors of estimated coefficients
    xi  <- rep(params[1], length(alpha))
    tau <- rep(params[2], length(alpha))
    
    result <- ifelse(abs(xi) < 0.0001, threshold - tau*log(alpha*n/k) ,
                                       threshold + tau/xi*((alpha*n/k)^-xi -1))
    return(result)

}

For GPD, there is a closed form for Expected Shortfall:

GPD_ES <- function(params, data, threshold, alpha = .05){
    # to vectorize the function we need vectors of estimated coefficients
    xi  <- rep(params[1], length(alpha))
    tau <- rep(params[2], length(alpha))    
    VAR <- GPD_VAR(params, data, threshold, alpha)
    
    result <- VAR/(1-xi)+ (tau-xi*threshold)/(1-xi)
    return(result)
}

The outputs are

GPD_VAR(c(0.15059, 0.86118), coredata(returns_xts), my_threshold, .01)
GPD_ES(c(0.15059, 0.86118), coredata(returns_xts), my_threshold, .01)

3.69068499489364
5.20816036412762

And we can check this with the output of fExtremes package function gpdRiskMeasures.

gpdRiskMeasures(fextreme_gpd_fit, prob= c(.99))

        p quantile shortfall
beta 0.99 3.690996 5.209194 

As you can see, our results are close to the output from fExtremes package.

Finally, a quick note on how to select the threshold u. This is more of an art than science. A common approach relies on Mean Excess Plots. The idea here is that for GPD, the mean excess value is a linear function of the threshold.

In fExtremes, there is a function mxfPlot that can do this for us.

mxfPlot(coredata(returns_xts), u = 0, cex = .5)

If you squint long enough, you might agree with me that a case can be made that mean excess is linear somewhere between 0.25-1.00. It seems like our initial choice of 0.85 wasn’t too bad. There are more methods that can assist us with the appropriate choice of threshold. You can consult some of the sources that are cited at the end of this post.

We covered a lot of the core concepts and introduced the main ideas behind EVT implementation. Hope you enjoyed.

Resources

• Good book with a dedicated chapter on EVT. I relied heavily on Chapter 6 for the discussion of MLE implementation in R. The author has an amazing blog that I discovered while finishing up this post. It is definitely worth checking out: https://freakonometrics.hypotheses.org/
• “An Introduction to Statistical Modeling of Extreme Values”
• Excellent article on EVT application in Finance.
• “Financial Risk Modelling and Portfolio Optimization in R” covers EVT in chapter 7. There is a lot of excellent material in this book in general, and is well worth the price.
• And just for fun, an interesting method to price catastrophe bonds based on EVT.

To highlight the benefits of building quantile regression models, I will contrast it with the ubiquitous linear regression model.

Regression models are used to study a relationship between a dependent (target) variable and some set of independent (explanatory) variables. We use explanatory variables to estimate what the mean/expected value of the target variables is conditioned on the observed explanatory variables. Here is the usual expression:

So what are we supposed to do when mean and variance are not enough. We need to get other information about the distribution. Higher order moments are important, e.g. skewness and kurtosis. Another way to describe the data is via its quantiles. For example, instead of the mean, we can use the median as the measure of central tendency. This often makes sense for skewed distributions. We can also examine interquartile range by looking at the difference between 75th and 25th percentile of the distribution.

This idea gave rise to the quantile regression model. Instead of modelling the mean/expected value of the target variable, let’s instead model some quantile of the data conditional on some explanatory variable. When we do that, we can estimate a set of models for some specified set of quantiles. Conditioning on explanatory variables, we will have a much more granular view of the distribution of the target variable. We can observe changes in the distribution (shape shift) at different values of the explanatory variable. This is incredibly useful when answering specific questions.

Let’s write an equation to solidify the concept. We have:

This is incredibly useful – we can analyze the impact of the explanatory variables on the full distribution of the target variable with fewer assumptions than in OLS.

Let’s use the example from quantreg vignette. See link here for more details.

require(quantreg)

data(engel)
attach(engel)
plot(income,foodexp,cex=.25,type="n",xlab="Household Income", ylab="Food Expenditure")
points(income,foodexp,cex=.5,col="blue")
abline(rq(foodexp~income,tau=.5),col="blue")
abline(lm(foodexp~income),lty=2,col="red") #the dreaded ols line
taus <- c(.05,.1,.25,.75,.90,.95)
for( i in 1:length(taus)){
    abline(rq(foodexp~income,tau=taus[i]),col="gray")
}

Here, we model household food expenditure as a function of income. The red line is the OLS regression line, the blue line is the median regression line, and the grey lines are the 5th, 10th, 25th, 75th, 90th, and 95th quantile lines. Notice that the quantiles capture shape shift of the distribution. At low income levels, the spread in the distribution is tight. For large incomes, household food expenditures have a very wide distribution. Also, notice that the difference is small between the 95th and 90th quantile, but the difference is big between the 5th and 10th quantile. There is so much more information here than usual OLS regression model.

I don’t want to do a deep dive of how the parameters are estimated. This is an interesting topic and more info is available in the references I provide at the end of the post. Let me just make two points. First, there is no closed form solution. We need our software to perform numerical routines and hope they converge to an optimal solution. Luckily this issue has been studied for a long time and reliable numerical routines have been implemented for us. Second, though it is not obvious, the solutions are derived by minimizing below function.

where

The idea is very similar to OLS where we estimate our parameters by minimizing some loss function. In the case of OLS we minimize sum of squared errors. In quantile regression, our cost function is asymmetric and depends on what quantile p we are estimating.

Now, time to get our hands dirty. Let’s work our way through an example that is very relevant in macro land these days. The Fed and markets are worried about the uncertainty about growth and market is pricing in insurance cuts from the Fed. Everyone seems to be in universal agreement that downside risks have increased this year and that it warrants a policy response from the Fed. So what we want to do here is estimate future GDP distribution conditioned on some economic variables and see if we can quantify the increase in uncertainty. There are some interesting papers that have been published by Fed researchers that tackle this problem directly. I provide the link at the end of the post, but click here if you want to have a look before proceeding.

The idea is that we can do the following:

1. Model the quantile of 1-year forward GDP as a function of current GDP and a broad measure of financial conditions. In the paper they use average GDP prints from current quarter out to 4 quarters ahead. I will instead model actual GDP print 4 quarters ahead. Also, the financials conditions index that we are using is released on a weekly basis. I use a 3-month average as an input into the model. More info about the index can be found here.
2. After estimating the quantiles, the authors fit a skewed Student-t distribution to get the full distribution. The authors solved for the 4 parameters of the distribution directly to match the 5th, 25th, 75th, and 95th quantile. In this post, we will estimate more quantile regression models and fit a Student-t by minimizing sum of squared residuals.
3. Once we have a full fitted distribution, we can calculate expected shortfall for GDP. We can then monitor the evolution of GDP tails risk conditional on current growth and financial conditions.

If you want to follow along, I have the data at the end of the blog. This data can be found online (publicly available) but I attached here for your convenience.

First, let’s fit the models for 5th, 10th, 15th, 25th, 35th, 50th, 65th, 75th, 85th, 90th, and 95th quantiles.

# fit quantile regression for 11 quantiles
my_taus <- c(.05, .10, .15, .25, .35, .50, .65, .75, .85, .90, .95)
my_qr_fit <- rq(GDP_Tp4 ~ GDP + NFCI, 
                tau = my_taus, 
                data = in_sample_data) # in_sampe_data is my raw data file

We can then look at the statistical significance of the estimated parameters. Please note that we need to run a bootstrap method to estimate the coefficients’ standard errors. Basically, we resample all the data with replacement many times (200 is the default). We estimate our coefficients using each sample. After we have bootstrapped estimates, we can calculate their standard errors

# check for statistical significance via bootstrap
my_summary <- summary(my_qr_fit, se = 'boot')
my_summary

Call: rq(formula = GDP_Tp4 ~ GDP + NFCI, tau = my_taus, data = in_sample_data)
 tau: [1] 0.05
 Coefficients:
             Value    Std. Error t value  Pr(>|t|)
 (Intercept) -2.34612  0.76494   -3.06705  0.00271
 GDP          0.32266  0.29841    1.08126  0.28190
 NFCI        -4.34642  0.81119   -5.35805  0.00000
 Call: rq(formula = GDP_Tp4 ~ GDP + NFCI, tau = my_taus, data = in_sample_data)
 tau: [1] 0.1
 Coefficients:
             Value    Std. Error t value  Pr(>|t|)
 (Intercept) -1.06657  0.80660   -1.32230  0.18876
 GDP         -0.03196  0.30036   -0.10642  0.91544
 NFCI        -4.45513  0.77950   -5.71534  0.00000
 Call: rq(formula = GDP_Tp4 ~ GDP + NFCI, tau = my_taus, data = in_sample_data)
 tau: [1] 0.15
 Coefficients:
             Value    Std. Error t value  Pr(>|t|)
 (Intercept) -0.77808  0.82633   -0.94161  0.34842
 GDP          0.03824  0.23594    0.16206  0.87155
 NFCI        -3.95458  0.89826   -4.40250  0.00002
 Call: rq(formula = GDP_Tp4 ~ GDP + NFCI, tau = my_taus, data = in_sample_data)
 tau: [1] 0.25
 Coefficients:
             Value    Std. Error t value  Pr(>|t|)
 (Intercept)  0.31032  0.80157    0.38713  0.69939
 GDP          0.02821  0.15774    0.17882  0.85840
 NFCI        -3.12373  1.13956   -2.74117  0.00713
 Call: rq(formula = GDP_Tp4 ~ GDP + NFCI, tau = my_taus, data = in_sample_data)
 tau: [1] 0.35
 Coefficients:
             Value    Std. Error t value  Pr(>|t|)
 (Intercept)  0.76880  0.75029    1.02468  0.30772
 GDP         -0.01704  0.15815   -0.10772  0.91441
 NFCI        -2.75392  0.90770   -3.03393  0.00300
 Call: rq(formula = GDP_Tp4 ~ GDP + NFCI, tau = my_taus, data = in_sample_data)
 tau: [1] 0.5
 Coefficients:
             Value    Std. Error t value  Pr(>|t|)
 (Intercept)  1.64539  0.65757    2.50221  0.01379
 GDP          0.02938  0.18618    0.15782  0.87488
 NFCI        -1.63139  0.84303   -1.93514  0.05549
 Call: rq(formula = GDP_Tp4 ~ GDP + NFCI, tau = my_taus, data = in_sample_data)
 tau: [1] 0.65
 Coefficients:
             Value    Std. Error t value  Pr(>|t|)
 (Intercept)  2.47940  0.56079    4.42127  0.00002
 GDP          0.09154  0.21980    0.41646  0.67787
 NFCI        -0.86369  0.82591   -1.04574  0.29793
 Call: rq(formula = GDP_Tp4 ~ GDP + NFCI, tau = my_taus, data = in_sample_data)
 tau: [1] 0.75
 Coefficients:
             Value    Std. Error t value  Pr(>|t|)
 (Intercept)  3.09624  0.39588    7.82120  0.00000
 GDP          0.14126  0.16218    0.87102  0.38561
 NFCI        -0.60402  0.62649   -0.96413  0.33706
 Call: rq(formula = GDP_Tp4 ~ GDP + NFCI, tau = my_taus, data = in_sample_data)
 tau: [1] 0.85
 Coefficients:
             Value    Std. Error t value  Pr(>|t|)
 (Intercept)  3.44024  0.30096   11.43101  0.00000
 GDP          0.15078  0.10782    1.39837  0.16476
 NFCI        -0.57208  0.43038   -1.32925  0.18647
 Call: rq(formula = GDP_Tp4 ~ GDP + NFCI, tau = my_taus, data = in_sample_data)
 tau: [1] 0.9
 Coefficients:
             Value    Std. Error t value  Pr(>|t|)
 (Intercept)  3.73520  0.24207   15.43012  0.00000
 GDP          0.22777  0.06743    3.37786  0.00101
 NFCI        -0.04564  0.41388   -0.11028  0.91238
 Call: rq(formula = GDP_Tp4 ~ GDP + NFCI, tau = my_taus, data = in_sample_data)
 tau: [1] 0.95
 Coefficients:
             Value    Std. Error t value  Pr(>|t|)
 (Intercept)  3.88112  0.16259   23.87047  0.00000
 GDP          0.22079  0.06693    3.29861  0.00130
 NFCI        -0.04759  0.39025   -0.12194  0.90317

There is a lot of output to digest. First, notice that for lower quantiles it is the financial conditions that are important, not current economic conditions. The alternative seems to be true for higher quantiles.

A much better way to view the output is via coefficient plots.

# coefficient plots
plot(my_summary)

From these plots we can immediately see that for lower quantiles, financial conditions index coefficients and their confidence bands are significantly different from zero. This is not the case for GDP coefficients.

Now that we have the fitted models, let’s extract the estimated quantiles.

# extract fitted quantiles based on model
my_predict <- predict(my_qr_fit)

We now need to be able to calibrate a Skewed Student-t (SST) distribution to these quantiles. To do this, we can set up a minimization routine. First, we will write an objective function that takes as arguments the 4 parameters of the SST distribution, the quantiles we are targeting, and the values of those quantiles. The function will return the sum of squared of the difference.

# objective function to fit
my_obj_fun <- function(my_dp, taus, target_quantile){
    # dp:                 trial a vector of length 4, elements represent location, scale, skew, and degrees of freedom
    # taus:               percentiles we are looking at (length of taus should match target_quantile)
    # target_quantile:    the empirical/model percentiles that we want GLD to match
    
    fitted_quantiles <- sn::qst(taus, dp = my_dp)
    not_inf <- !is.infinite(fitted_quantiles)
    result <- sqrt(sum((target_quantile[not_inf] - fitted_quantiles[not_inf])^2))*100
    return(result)
}

We will then write a function that uses R’s optimization routine to minimize our objective function:

# minimization routine
my_solver <- function(target_quantile, init_values, taus){
    result <- optim(init_values, 
                 my_obj_fun, 
                 taus = taus,  
                 target_quantile = target_quantile, 
                 method = 'Nelder-Mead', 
                 control = list(maxit = 100, abstol = .000001))$par
    
    return(result)
}

Now we will first fit an SST distribution to all our GDP observations to get a set of good starting values to feed into our minimization routine.

require(sn)
SST_Fit <- sn::st.mple(y=in_sample_data[,'GDP'])
initial_values <- SST_Fit$dp

And let’s plot the fitted distribution along with the histogram of the actual distribution of quarterly GDP prints.

hist(in_sample_data[,'GDP'],
 col='peachpuff', 
 border='black',
 prob = TRUE, # show densities instead of frequencies
 xlab = 'GDP',
 xlim=c(-5,6),
 breaks=20,
 main = '')

lines(x = seq(-5,6, by = .05), 
      y = dst(seq(-5,6, by = .05), dp=initial_values),
 lwd = 3, # thickness of line
 col = 'steelblue')

Let’s try this now to calibrate a SST distribution to the first set of fitted quantiles from our regression model.

# try solver on one observation
test_param <- my_solver(target_quantile = my_predict[1,], 
                        initial_values, 
                        taus = my_taus)

And let’s see quantiles from our calibrated distribution based on the fitted parameters.

sn::qst(my_taus, dp = test_param)

-1.82167377991434  -0.897478033842928  -0.317201996236203
0.488771386216558  1.09817828817344  1.87586807146716  2.62893251134214
3.18728127939248  3.88910109003964  4.3725835716495  5.11220664005463

my_predict[1,]

-1.08828972132954  -1.1917950846466  -0.629463584621562
0.41992958066738  0.702009508800799  1.75977359546504  2.83628037426853
3.64709619438839  4.02819298775264  4.62350035122945  4.74217692823193

And one more time for another set of estimated quantiles via quantile regression:

# try solver on one observation
test_param <- my_solver(target_quantile = my_predict[10,], 
                        initial_values, 
                        taus = my_taus)

sn::qst(my_taus, dp = test_param)

0.361081016667099  0.976569681870633  1.35985579676273  1.88170447048113
2.26245622070612  2.72305723000699  3.13679995549024  3.42278780682794
3.75976486690277  3.97940161356939  4.29917698946775

my_predict[10,]

0.297148107584487  1.20757328502942  1.32048689850673  1.96558089844999
2.17781950506346  2.52744469760649  3.03755265206431  3.57928464901055
3.91811899382793  4.03221608883146  4.17076043833115

Now, let’s loop through all our observations and fit an SST distribution at each point in time.

fitted_sst_params <- t(apply(my_predict, 
                             MARGIN = 1, 
                             FUN = my_solver, 
                             init_values = initial_values, 
                             taus = my_taus))

We can now compute GDP at 5% VaR. This is basically the 5th quantile that can be calculated as below.

N <- nrow(my_predict)
GDP_VAR <- rep(NA, N)
for (ii in 1:N){
    GDP_VAR[ii] <- sn::qst(p = .05, dp = fitted_sst_params[ii,])
}

Now, let’s plot a time series of GDP at 5% VaR. In our plot, the 5th quantile from the fitted SST is plotted as a line. I included grey dots that come straight from the fitted quantile regression for reference also. Any difference between the two is simply a function of our fitting error earlier.

plot(in_sample_data[,'Date'], 
     GDP_VAR, 
     type = 'l', 
     xlab='Date', 
     ylab=' GDP VaR')

points(in_sample_data[,'Date'],my_predict[,1], pch = 19, col = rgb(.25,.25,.25,.25))

Next, calculate the expected shortfall.

# function for Expected Shortfall (Conditional VAR) (numerically integrate from 0 to alpha)
SST_ES <- function(params, alpha = .05){
    result <- 1/(alpha) * integrate(sn::qst, lower = .00001, upper = alpha, stop.on.error = TURE, dp = params)$value
    return(result)
}

# calculate ES for each quarter
GDP_ES <- apply(fitted_sst_params, MARGIN = 1, FUN = SST_ES)

# Plot ES on a rolling-base
plot(in_sample_data[,'Date'], GDP_ES, type = 'l', xlab='Date', ylab=' GDP ES', ylim = c(-20,5))

Notice that we get close to the results of Adrian et al. (2019).

Note that our results are not identical because we have defined the problem slightly differently. Here we are modelling the distribution of the quarterly GDP print four quarters in the future, while Adrian et al. (2019) are modelling the distribution of the average of four furture quarterly GDP readings.

Now, I want to demonstrate a 3-D plot of how GDP distribution evolves over time according to our model.

# data for plot
GDP_Plot <- seq(-25, 25, by = .5)
Prob_Plot <- matrix(NA, nrow = length(in_sample_data[,'Date']), ncol = length(GDP_Plot))

for (ii in 1:length(in_sample_data[,'Date'])){
    Prob_Plot[ii,] <- sn::dst(GDP_Plot, dp = fitted_sst_params[ii,])
}


require(grDevices)
persp(x = 1:115,
      y = GDP_Plot,
      z = Prob_Plot,
      xlab = 'Date',
      ylab = 'GDP t+4q',
      zlab = 'Prob',
      xlim = c(1,115),
      ylim = c(-25,25),
      zlim = c(0,1.1),
      theta = 75,
      phi = 15,
      col = 'grey9',
      #col =  color[facetcol],
      border='grey40',
      expand = 0.5)

Finally, we can also generate scenarios using our model. We first create a data frame to hold our scenarios, then feed them into our quantile regression to get the estimated quantiles, and finally, we fit the skewed student t distribution to the quantiles.

# scenario data
scenario_data <- data.frame(GDP = c(-1.0, 1.5, 3.0),
                            NFCI = c(1.5, 0.0, -.50),
                            stringsAsFactors = FALSE)

#fitted quantiles based on model
scenario_quantiles <- predict(my_qr_fit, newdata = scenario_data)

#skewed student t parameters for scenarios
scenario_sst_params <- t(apply(scenario_quantiles, 
                             MARGIN = 1, 
                             FUN = my_solver, 
                             init_values = initial_values, 
                             taus = my_taus))

Now we can plot the distributions.

# plot 
gdp_grid <- seq(-20,15, by = .25)
plot(gdp_grid, sn::dst(gdp_grid, dp = scenario_sst_params[1,]), type= 'l', lty = 1,
     xlab='Prob', ylab= 'GDP t+4q',  ylim=c(0, .4))
polygon(gdp_grid, sn::dst(gdp_grid, dp = scenario_sst_params[1,]),col=rgb(.5, .95, .75, .15))

lines(gdp_grid, sn::dst(gdp_grid, dp = scenario_sst_params[2,]))
polygon(gdp_grid, sn::dst(gdp_grid, dp = scenario_sst_params[2,]),col=rgb(.7, .1, .4, 0.15))


lines(gdp_grid, sn::dst(gdp_grid, dp = scenario_sst_params[3,]))
polygon(gdp_grid, sn::dst(gdp_grid, dp = scenario_sst_params[3,]),col=rgb(1, 0, 1, .15))

legend("topleft", legend=c('Scn1', 'Scn2', 'Scn3'),
       col=c('seagreen3', 'violetred4', 'violet'), lty = 1, inset=.02, cex = .6, y.intersp=2.5)

I hope I convinced you that quantile regression is a powerful tool. It can be very effective for capturing conditional distribution of some variable of interest. Good luck!

Resources

• Adrian, Tobias, et al. Vulnerable Growth – FEDERAL RESERVE BANK of NEW YORK, www.newyorkfed.org/research/staff_reports/sr794
• Loria, et al. “Assessing Macroeconomic Tail Risk.” SSRN, 19 Apr. 2019, papers.ssrn.com/sol3/papers.cfm?abstract_id=3375006
• FINANCIAL CONDITIONS AND GROWTH AT RISK – IMF. www.imf.org/en/Publications/GFSR/Issues/2017/09/27/~/media/Files/Publications/GFSR/2017/October/chapter-3/Documents/c3.ashx?la=en
• “National Financial Conditions Index (NFCI).” National Financial Conditions Index (NFCI) – Federal Reserve Bank of Chicago, www.chicagofed.org/publications/nfci/index
• Koenker, Roger. “Quantile Regression [R Package Quantreg Version 5.41].” The Comprehensive R Archive Network, Comprehensive R Archive Network (CRAN), cran.r-project.org/web/packages/quantreg/index.html

# data for the model
 in_sample_data <- read.csv(header = TRUE, 
                             text = 'Date, GDP, NFCI, GDP_Tp4 
                                     30-09-1989, 3.9, 0.000127008, 1.7
                                     31-12-1989, 2.7, -0.036629133, 0.6
                                     31-03-1990, 2.8, -0.090870083, -1
                                     30-06-1990, 2.4, -0.182160633, -0.5
                                     30-09-1990, 1.7, -0.02811965, -0.1
                                     31-12-1990, 0.6, 0.277148608, 1.2
                                     31-03-1991, -1, 0.091021233, 2.9
                                     30-06-1991, -0.5, -0.288714925, 3.2
                                     30-09-1991, -0.1, -0.480414508, 3.7
                                     31-12-1991, 1.2, -0.51906455, 4.4
                                     31-03-1992, 2.9, -0.564142033, 3.3
                                     30-06-1992, 3.2, -0.718176483, 2.8
                                     30-09-1992, 3.7, -0.685154442, 2.3
                                     31-12-1992, 4.4, -0.59803585, 2.6
                                     31-03-1993, 3.3, -0.755371775, 3.4
                                     30-06-1993, 2.8, -0.799939058, 4.2
                                     30-09-1993, 2.3, -0.896682492, 4.3
                                     31-12-1993, 2.6, -0.824130317, 4.1
                                     31-03-1994, 3.4, -0.821335817, 3.5
                                     30-06-1994, 4.2, -0.685582642, 2.4
                                     30-09-1994, 4.3, -0.727856025, 2.7
                                     31-12-1994, 4.1, -0.567921583, 2.2
                                     31-03-1995, 3.5, -0.608324375, 2.6
                                     30-06-1995, 2.4, -0.671981133, 4
                                     30-09-1995, 2.7, -0.648721525, 4.1
                                     31-12-1995, 2.2, -0.696705942, 4.4
                                     31-03-1996, 2.6, -0.681872408, 4.3
                                     30-06-1996, 4, -0.649720717, 4.3
                                     30-09-1996, 4.1, -0.677856925, 4.7
                                     31-12-1996, 4.4, -0.692815608, 4.5
                                     31-03-1997, 4.3, -0.661724692, 4.9
                                     30-06-1997, 4.3, -0.656140467, 4.1
                                     30-09-1997, 4.7, -0.655591383, 4.1
                                     31-12-1997, 4.5, -0.532616742, 4.9
                                     31-03-1998, 4.9, -0.6195233, 4.8
                                     30-06-1998, 4.1, -0.621981908, 4.7
                                     30-09-1998, 4.1, -0.442741983, 4.7
                                     31-12-1998, 4.9, -0.136276067, 4.8
                                     31-03-1999, 4.8, -0.345186308, 4.2
                                     30-06-1999, 4.7, -0.429547225, 5.3
                                     30-09-1999, 4.7, -0.2570176, 4.1
                                     31-12-1999, 4.8, -0.277550325, 3
                                     31-03-2000, 4.2, -0.293507592, 2.3
                                     30-06-2000, 5.3, -0.176782192, 1.1
                                     30-09-2000, 4.1, -0.283548017, 0.5
                                     31-12-2000, 3, -0.220909875, 0.2
                                     31-03-2001, 2.3, -0.257047708, 1.3
                                     30-06-2001, 1.1, -0.3745106, 1.3
                                     30-09-2001, 0.5, -0.3794091, 2.2
                                     31-12-2001, 0.2, -0.275062117, 2.1
                                     31-03-2002, 1.3, -0.384926792, 1.8
                                     30-06-2002, 1.3, -0.518591775, 2
                                     30-09-2002, 2.2, -0.28980035, 3.3
                                     31-12-2002, 2.1, -0.305308933, 4.3
                                     31-03-2003, 1.8, -0.390609775, 4.3
                                     30-06-2003, 2, -0.592714233, 4.2
                                     30-09-2003, 3.3, -0.579556617, 3.4
                                     31-12-2003, 4.3, -0.636004083, 3.3
                                     31-03-2004, 4.3, -0.728082242, 3.9
                                     30-06-2004, 4.2, -0.70020205, 3.6
                                     30-09-2004, 3.4, -0.676979208, 3.5
                                     31-12-2004, 3.3, -0.708622083, 3.1
                                     31-03-2005, 3.9, -0.700670808, 3.4
                                     30-06-2005, 3.6, -0.618226142, 3.1
                                     30-09-2005, 3.5, -0.624520042, 2.4
                                     31-12-2005, 3.1, -0.621440575, 2.6
                                     31-03-2006, 3.4, -0.655807208, 1.5
                                     30-06-2006, 3.1, -0.665672808, 1.8
                                     30-09-2006, 2.4, -0.65326455, 2.2
                                     31-12-2006, 2.6, -0.684736175, 2
                                     31-03-2007, 1.5, -0.70984855, 1.1
                                     30-06-2007, 1.8, -0.659238383, 1.1
                                     30-09-2007, 2.2, -0.175482867, 0
                                     31-12-2007, 2, 0.095068958, -2.8
                                     31-03-2008, 1.1, 0.493424692, -3.3
                                     30-06-2008, 1.1, 0.439167758, -3.9
                                     30-09-2008, 0, 0.722941642, -3
                                     31-12-2008, -2.8, 2.342395083, 0.2
                                     31-03-2009, -3.3, 1.827469433, 1.7
                                     30-06-2009, -3.9, 1.027353242, 2.8
                                     30-09-2009, -3, 0.394184042, 3.2
                                     31-12-2009, 0.2, 0.118619033, 2.6
                                     31-03-2010, 1.7, -0.125451367, 1.9
                                     30-06-2010, 2.8, -0.095603975, 1.7
                                     30-09-2010, 3.2, -0.114889233, 0.9
                                     31-12-2010, 2.6, -0.25043795, 1.6
                                     31-03-2011, 1.9, -0.341584208, 2.7
                                     30-06-2011, 1.7, -0.358413492, 2.4
                                     30-09-2011, 0.9, -0.007418492, 2.5
                                     31-12-2011, 1.6, 0.081366008, 1.5
                                     31-03-2012, 2.7, -0.227845667, 1.6
                                     30-06-2012, 2.4, -0.232090867, 1.3
                                     30-09-2012, 2.5, -0.341464892, 1.9
                                     31-12-2012, 1.5, -0.480014042, 2.6
                                     31-03-2013, 1.6, -0.562598125, 1.5
                                     30-06-2013, 1.3, -0.605581217, 2.6
                                     30-09-2013, 1.9, -0.592675233, 3
                                     31-12-2013, 2.6, -0.717339183, 2.7
                                     31-03-2014, 1.5, -0.7624705, 3.8
                                     30-06-2014, 2.6, -0.790339467, 3.4
                                     30-09-2014, 3, -0.774638708, 2.4
                                     31-12-2014, 2.7, -0.71241615, 2
                                     31-03-2015, 3.8, -0.625803008, 1.6
                                     30-06-2015, 3.4, -0.698488767, 1.3
                                     30-09-2015, 2.4, -0.622464742, 1.5
                                     31-12-2015, 2, -0.600945333, 1.9
                                     31-03-2016, 1.6, -0.501130675, 1.9
                                     30-06-2016, 1.3, -0.561193167, 2.1
                                     30-09-2016, 1.5, -0.565286533, 2.3
                                     31-12-2016, 1.9, -0.625796542, 2.5
                                     31-03-2017, 1.9, -0.700084275, 2.6
                                     30-06-2017, 2.1, -0.738209208, 2.9
                                     30-09-2017, 2.3, -0.788917608, 3
                                     31-12-2017, 2.5, -0.83847425, 3
                                     31-03-2018, 2.6, -0.781312675, 3.2',
                             stringsAsFactors = FALSE)

To begin, I will introduce the Generalized Lambda Distribution (GLD).  There is ample evidence to show that this is an excellent contender for modelling various financial returns.  At the end  of the post I reference two papers that test GLD on FX and Equity returns.

GLD is defined in terms of the inverse cumulative distribution function (cdf):

The only restrictions we have on the parameters is that has to be positive. (This is not unusual – for example, variance in the normal distribution has to be greater than zero).

In R, there is a great package GLDEX that has all the functionality that we need to calibrate the FMKL GLD distribution to data. We will work with SPX data from 2004-2019 and calibrate the distribution parameters to the full data set.

require(quantmod)
require(GLDEX)

SPY_prices <- Ad(getSymbols("SPY", from = "2004-01-01", to = "2019-01-01",
                 auto.assign = FALSE))
returns_xts <- diff(log(SPY_prices), k =1)[-1,]*100
full_sample_fit <- fun.RMFMKL.ml(data = returns_xts)  # MLE fit

fun.plot.fit(full_sample_fit, returns_xts, xlab = 'SPY Return',
             name = 'FMKL.ML', nclass = 75, param.vec = 'fmkl')

We can also check the qq plot to see how well the distribution fits the data.

qqplot.gld(returns_xts, full_sample_fit, 'fmkl', name = 'RMFMKL.ML')

The nice thing is that GLD is already specified as an inverse cumulative normal distribution so we can very easily use it to check Value at Risk (VaR) at a specified level .  For example, VaR at = 0.01 is -3.8%. The 99th percentile of the distribution is 3.06%.

qgl(.01, lambda1= full_sample_fit)
qgl(.99, lambda1= full_sample_fit)

-3.77901522473645
3.05837784742503

One thing the package does not give us is the expected shortfall formula (ES).  ES is the conditional VaR. i.e. the expected loss given that our VaR has been exceeded.

We can implement the function below.

ES_FMKL_GLD <- function(alpha, lambda1, lambda2=NULL, lambda3=NULL, lambda4=NULL){
    if (length(lambda1) > 1) {
        lambda2 <- lambda1[2]
        lambda3 <- lambda1[3]
        lambda4 <- lambda1[4]
        lambda1 <- lambda1[1]
    }
    lambda1 + 1/(alpha*lambda2*lambda3) * ((alpha^(lambda3+1))/(lambda3+1)-alpha) - 
              1/(alpha*lambda2*lambda3) * ((1-(1-alpha)^(lambda4+1))/(lambda4+1) - alpha)
}

Let’s check it against numerical integration:

alpha <- .01
round(ES_FMKL_GLD(alpha, lambda1 = full_sample_fit),4)
round(1/(alpha) * integrate(GLDEX::qgl, lower = 0, upper = alpha, 
                           stop.on.error = FALSE, lambda1 = full_sample_fit)$value,4)

-6.4877
-6.4873

They seem close enough.

Now let’s check what the 5% VaR is and the expected shortfall. Notice that while the VaR is not that large (-1.66%), the expected shortfall is much more significant (-3.18%). This is a reason why VaR is not a great measure of risk.

qgl(.05, lambda1= full_sample_fit)
ES_FMKL_GLD(.05, lambda1 = full_sample_fit)

-1.65945064222568
-3.18137496967441

For the next exercise we will set up an out-of-sample procedure where we use a rolling 2-year window to calibrate the GLD parameters, then we estimate 0.1%, 1%, 5%, 95%, 99%, 99.9% VaR for the next day, and finally we compare it to subsequent day’s return. Based on all the observations we would expect 5% and 95% VaR to be exceeded in approximately 5% of the sample.

First, let’s set up the data frame that will hold our results.

N <- length(returns_xts)
lookback <- 251*2
sample_end <- lookback:N
sample_start <- 1:length(sample_end)
alpha1 <- .95
alpha2 <- .99
alpha3 <- .999

output_df <- data.frame(Date = as.Date(rep(NA, N-lookback)),
                       Lambda1 = numeric(N-lookback),
                       Lambda2 = numeric(N-lookback),
                       Lambda3 = numeric(N-lookback),
                       Lambda4 = numeric(N-lookback),
                       VaR_001 =  numeric(N-lookback),
                       VaR_01 =  numeric(N-lookback),
                       VaR_05 =  numeric(N-lookback),
                       VaR_95 =  numeric(N-lookback),
                       VaR_99 =  numeric(N-lookback),
                       VaR_999 =  numeric(N-lookback),
                       Rtrn_Tp1 =  numeric(N-lookback),
                       stringsAsFactors = FALSE)

Below is the main loop for the backtest. Please note that here I specify VaR 95 at the 5th percentile of returns because it is typical in finance to think of VaR in terms of losses. When one discusses 95% VaR it is the 5th percentile of the distribution that we are interested in. In the loop I use (1-alpha) as an argument to the quantile function.

# Main Loop
for (ii in 1:(N-lookback)){
    subset_sample <- coredata(returns_xts[sample_start[ii]:sample_end[ii],]) # extract the window
    my_gld_fit <- fun.RMFMKL.ml(data = subset_sample)  # calibrate parameters
    my_VaR <- qgl(c(c(alpha3, alpha2, alpha1), 1- c(alpha1, alpha2, alpha3)), lambda1 = my_gld_fit)  # calculate VaR
    
    # extract results
    output_df[ii, 'Date'] <- index(returns_xts[sample_end[ii]]) 
    output_df[ii, c(2,3,4,5)] <- my_gld_fit
    output_df[ii, c(6,7,8,9,10,11)] <- my_VaR
    output_df[ii, 'Rtrn_Tp1'] <- returns_xts[sample_end[ii]+1,1]
}

tail(output_df)

As a very quick check we simply count the number of times the observed t+1 return exceeded the estimated VaR. You can see that the fitted distribution does a decent job of matching our expectations. 5% VaR was exceeded 5.25% of the time. 95% VaR was exceeded 5.78% of the time.

round(sum(output_df[,'Rtrn_Tp1'] > output_df[,'VaR_001'])/length(returns_xts)*100,2)
round(sum(output_df[,'Rtrn_Tp1'] > output_df[,'VaR_01'])/length(returns_xts)*100,2)
round(sum(output_df[,'Rtrn_Tp1'] > output_df[,'VaR_05'])/length(returns_xts)*100,2)

round(sum(output_df[,'Rtrn_Tp1'] < output_df[,'VaR_95'])/length(returns_xts)*100,2)
round(sum(output_df[,'Rtrn_Tp1'] < output_df[,'VaR_99'])/length(returns_xts)*100,2)
round(sum(output_df[,'Rtrn_Tp1'] < output_df[,'VaR_999'])/length(returns_xts)*100,2)

0.21
1.35
5.25
5.78
1.38
0.11

Below I plotted daily returns and the 99% VaR as a red line.

plot(output_df[,'Date'], 
     output_df[,'Rtrn_Tp1'], 
     cex = .5, 
     pch = 19, 
     col = rgb(red = .275, green = .51, blue = .706, alpha =.25),
     ylab = 't+1 Return',
     xlab = 'Date')

lines(output_df[,'Date'], output_df[,'VaR_99'], col = 'red', lwd = 2)

This was a quick guide on how to implement a useful flexible distribution to model asset returns. In future posts I will introduce and explore other flexible distributions. Stay tuned.

Useful Resources:
1) Corlu, Canan G., and Alper Corlu. “Modelling Exchange Rate Returns: Which Flexible Distribution to Use?” Quantitative Finance, vol. 15, no. 11, 2014, pp. 1851–1864., doi:10.1080/14697688.2014.942231. Link

2) Corlu, Canan Gunes, et al. “Empirical Distributions of Stock Returns: A Comparison.” SSRN Electronic Journal, 2015, doi:10.2139/ssrn.2622672. Link

The Algorithm:

Since Nelder-Mead (NM) is a derivatives-free algorithm it can be applied to optimization problems for which the objective function f (x) is not differentiable. As long as we can evaluate f (x) at any point x we can apply the algorithm.  The downside is that NM is a heuristic technique which means that it uses a practical search method that is not guaranteed to converge to the optimal solution.

In this section I will discuss the steps of the algorithm.  Let’s define the problem as:

This simply means that we are trying to minimize some multivariate function f () that has n dimensional input vector x and where each variable in x is continuous.  I will later discuss ways in which we can modify the algorithm to handle problems where some variables in x are discrete or are bounded.

A quick side note, NM is often referred to as Nelder-Mead Simplex, Downhill Simplex, or Amoeba algorithm.  This is because at any point in time we are working with n+1 points (ie one point more than the dimension of the problem).  Each point that is joined by a line (or edge) is called a vertex.  When we plot these points they form a simplex.  In one dimension a simplex is a line, in two dimensions a triangle, and in three dimensions it is a tetrahedron.

I will use an example in two dimensional space so that we can visually inspect the geometry of the algorithm.  We begin by selecting some n+1 points (three points in this example) to form an initial simplex (triangle). Our vertices are(remember that each x is of dimension n):

We sort these vertices in increasing order of the objective function values:

We then compute a centroid which is the mean over all the vertices excluding the worst (omit vertex that has the highest objective function value).  For our example I will use the following data:

I do not specify an objective function but since we assume the data has been ordered f (x_1) achieves the minimum at this stage.

Our centroid can be calculated as:

After we have our centroid we compute a reflection of the worst vertex:

Usually rho is set to 1, so we get:

Below is a plot of our initial simplex in grey and the new simplex formed by our reflection point (highlighted in blue).

Label the new point x_r.  We now evaluate our objective function at this new point.  If the objective function at point x_r, f(x_r), is lower than our current minimum at x_1 then we compute an expansion.  This just means we continue in the same direction to see if we can improve on x_r.

Expansion is given as:

Chi is usually set to 2, using that we have:

Below is a plot of the new simplex.  Expansion is highlighted in red.  I kept the reflection point so you can see the geometry.

We can now evaluate our objective function at point x_e.  If we further improved on f(x_r),  ie the objective function at x_e is lower than at x_r, then we remove the worst vertex x_n+1 and replace it with x_e.  If f(x_e) is not an improvement on f (x_r) then we replace the worst vertex with x_r.

So far we have taken care of the case where at the very least our reflection operation improved our solution.  Now let’s consider what can be done if f(x_r) was not an improvement on the current minimum objective function value f(x_1).  We can have a situation where f( x_r) is worse than f(x_1) but is better than the second worst f(x_n).  In this case we still improved on the worst vertex so we would accept x_r and replace x_n+1 with it.

Now what if f(x_r) is an improvement on the worst vertex but is not better than the second worse vertex.  In this case NM considers a contraction outward.

Out-Contraction is given as:

Parameter gamma is usually assigned a value of .5, so we end up with:

We are basically dampening the reflection operation to see if we moved too far initially.

If our new objective function value at point x_o is less than f(x_r) then we accept this change and replace the worst vertex x_n+1 with x_0.  Below is a plot of the contraction in red and can be compared to the reflection point.

If out contraction did not improve on reflection then we will shrink the simplex.  Shrinking the simplex means we move all vertices closer to the best vertex x_1.  Shrink formula is below:

Sigma is typically set equal to .5.  Therefore, for x_2 we get:

And for x_3:

Below is a plot of our original simplex and the new shrunk points that are joined by a purple edge.

In other cases we consider a contraction inward.  The formula is given below:

You can see that in-contraction is a move of the worst vertex toward the centroid.  Below plot includes out-contraction in red and in-contraction in green.

Finally, if at this stage we have not improved on the worst vertex we again shrink the simplex.

At this point I probably have lost most of you.  For those who are still with me I created a flowchart of the algo below.  Its a busy chart but I think its worthwhile to spend some time inspecting it.  The blue rectangles are points in the algorithm where we have to evaluate the objective function.  The green circles is where we swap out the worst vertex for some improvement. Finally, the red rectangles are our test conditions where we search for an improvement of the initial simplex.

Modifications:

The NM algorithm has a few drawbacks besides the fact that it is a heuristic approach.  The first issue is that the input vector is continuous.  In many cases we want to search in a bounded space.  For an example in finance, some SABR or Heston model parameters are bounded, therefore we don’t want to search all the parameter space when we are calibrating the models.  In Heston and SABR we have parameter rho as an input.  Rho captures spot-vol correlation in the model and should be in range (-1,1).  An example of search in a bounded space in the field of machine learning is when we are using NM to tune a hyper-parameter of our model.  For example, when calibrating a support vector machine we may want to find an optimal Cost parameter that minimizes some cross validation error.  Our tuning parameter must be greater than zero.  So C should be in range (0,inf).

One way to deal with this is by filtering our objective function.  We can assign a large value to f(x) when any x is outside of our bounded search space.  So in the case of SABR or Heston, when NM calculates a reflection and takes the value of Rho outside of the bounded space we set the value of the objective function to some high value.  This will force the algorithm to perform a contraction and keep us in the search space we desire.  This may work but will increase the number of iterations it takes for the algorithm to converge, especially if the minimum is close to the boundary value of Rho.

A more preferred way to deal with the issue is to map our constrained optimization problem to an unconstrained one.  For example, If Rho needs to be in range (-1,1) we can amend our objective function to be a function of some variable x where:

So our minimization problem becomes:

We can simply take the inverse of the tanh function to get the value of x.

For boundaries that deal with range (0,inf), like our Cost parameter in case of SVM  we can use:

To get back our x value we can take the natural log of Cost.

Another problem that we can often encounter is that our variable x is discrete not continuous.  Again let’s take an example from machine learning.  Some hyper-parameters are discrete.  For instance, the number of nodes or hidden layers in a neural network must be discrete.  One possible solution to this problem is to round x to the nearest integer value after every change we make to the simplex.

There are other issues with NM such as collapse of a simplex, premature convergence on local (not global) minimum, failure to converge (stuck in a plateau).  There have been many modifications suggested to deal with these issues.  Thing like multiple restarts with different starting simplexes and approximation to gradients around the simplex to infer more information about the objective function have been tried and work well in many fields.

VBA implementation:

Since I cannot upload vba enabled workbooks I am attaching a macro-free workbook and the code is at the end of the post.

1. Open Excel workbook hit Alt+F11 to open VB Editor.  Right click on ThisWorkbook and Insert Module.  NelderMead_DownhillSimplex
2. Highlight, Copy & Paste the code at the end of this post in its entirety into the Module.

In the workbook I have an example of a where we try to minimize the banana function using starting values (-1.2,1). Formula Nelder_Mead_Min is in range B7:D9.  This is an array formula so you need to hit Ctrl_Shift_Enter.  After 76 iterations we converge on an answer that is close to the true min of (1,1)

Useful Resources:

• VBA Developer’s Handbook is an excellent book.  It was the source for the bubble sort algo that I modified http://ca.wiley.com/WileyCDA/WileyTitle/productCd-078215333X.html
• Wikipedia entry for Bubble Sort Algo https://en.wikipedia.org/wiki/Bubble_sort
• Rosenbrock Function: https://en.wikipedia.org/wiki/Rosenbrock_function#An_example_of_optimization
• Interesting list of optimization test functions http://www.sfu.ca/~ssurjano/optimization.html  and https://arxiv.org/pdf/1308.4008.pdf
• Scilab (Open source numerical computation software) has an excellent piece on Nelder-Mead “Nelder-Mead User’s Manual” by Michael Baudin https://wiki.scilab.org/The%20Nelder-Mead%20Component?action=AttachFile&do=view&target=neldermead.pdf
• Excellent book “Option Pricing Models and Volatility” that includes VBA implementation of Heston model that uses Nelder Mead for calibration http://ca.wiley.com/WileyCDA/WileyTitle/productCd-0471794643.html
• Christopher Brooks has a short and interesting paper on the NM algo https://www.cs.usfca.edu/~brooks/papers/amoeba.pdf

Copy Code Starting Below:

Option Explicit
Option Base 1

Function ObjectivFun_Banana(ArgIn)
    Dim x1 As Double:   x1 = ArgIn(1)
    Dim x2 As Double:   x2 = ArgIn(2)
    
    ObjectivFun_Banana = (1 - x1) ^ 2 + 100 * (x2 - x1 ^ 2) ^ 2
End Function

Public Function Nelder_Mead_Min(function_name As String, function_parameters As Variant, Optional MaxIterations As Integer = 150)

'Factors for reflection, expansion, contraction, and shrinkage
    Dim rho As Double:      rho = 1
    Dim chi As Double:      chi = 2
    Dim gamma As Double:    gamma = 0.5
    Dim sigma As Double:    sigma = 0.5

'xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx’
'x Variables that will hold function values for
'x reflection, expansion, and both contraction points
'xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx’
    Dim ReflectionValue As Double:          Dim ExpansionValue As Double
    Dim OutContractionValue As Double:      Dim InContractionValue As Double
    Dim ShrinkFlag As Boolean

    Dim ParamCount As Integer:  ParamCount = CInt(function_parameters.Cells.Count) 'dimension of Fun

'xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx’
'x Arrays that will hold new vertices
'xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx’
    Dim x_bar() As Double:          ReDim x_bar(ParamCount)
    Dim x_r() As Double:            ReDim x_r(ParamCount)
    Dim x_e() As Double:            ReDim x_e(ParamCount)
    Dim x_o() As Double:            ReDim x_o(ParamCount)
    Dim x_i() As Double:            ReDim x_i(ParamCount)

    Dim ii As Integer: ii = 1 ' counter
    Dim jj As Integer: jj = 1 ' counter

    Dim Temp_x_vec() As Double: ReDim Temp_x_vec(ParamCount)

'xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx’
'x  creates a  ParamCount+1 x ParamCount+1 matrix called Simplex_Matrix
'x  1st column is function value  while remaining columns are the vertecise
'x  page 20 https://wiki.scilab.org/The%20Nelder
'x  -Mead%20Component?action=AttachFile&do=view&target=neldermead.pdf
'xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
    Dim Simplex_Matrix() As Double:  ReDim Simplex_Matrix(ParamCount + 1, ParamCount + 1)
    Dim delta_u As Double: delta_u = 0.05
    Dim delta_z As Double: delta_z = 0.0075
    
    For ii = 1 To (ParamCount + 1)
        For jj = 1 To ParamCount
            If ii = 1 Then ' first row uses supplied starting values
                Simplex_Matrix(ii, jj + 1) = function_parameters(jj)
            ElseIf (jj = ii - 1) And (function_parameters(jj) <> 0) Then
                Simplex_Matrix(ii, jj + 1) = function_parameters(jj) * (1 + delta_u)
            ElseIf (jj = ii - 1) And (function_parameters(jj) = 0) Then
                Simplex_Matrix(ii, jj + 1) = delta_z
            Else
                Simplex_Matrix(ii, jj + 1) = function_parameters(jj)
            End If
            Temp_x_vec(jj) = Simplex_Matrix(ii, jj + 1) ' temp vector of x to evaluate ObjFun
        Next jj
        Simplex_Matrix(ii, 1) = Run(function_name, Temp_x_vec) ' eval ObjFun
    Next ii

' sort simplex_matrix based on FunObj
    Simplex_Matrix = BubbleSort(Simplex_Matrix)

    Dim TotalIt As Integer: TotalIt = 1 ' iteration count
    Dim Epsilon As Double:  Epsilon = 0.000001 ' convergence criteria
' main loop
    Do While Abs(Simplex_Matrix(1, 1) - Simplex_Matrix(ParamCount + 1, 1)) > Epsilon And TotalIt < MaxIterations

' compute centroid
        For jj = 1 To ParamCount
            x_bar(jj) = 0
            For ii = 1 To ParamCount ' sum all x excluding best
                x_bar(jj) = x_bar(jj) + Simplex_Matrix(ii, jj + 1)
            Next ii
            x_bar(jj) = x_bar(jj) / ParamCount ' compute average
        Next jj

        For jj = 1 To ParamCount ' compute reflection x_r
            x_r(jj) = (1 + rho) * x_bar(jj) - rho * Simplex_Matrix(ParamCount + 1, jj + 1)
        Next jj
        ReflectionValue = Run(function_name, x_r) ' compute f(x_r)

        If ReflectionValue < Simplex_Matrix(1, 1) Then ' if f(x_r)<f(x_1) then
            For jj = 1 To ParamCount ' compute expansion x_e
                x_e(jj) = (1 + rho * chi) * x_bar(jj) - rho * chi * Simplex_Matrix(ParamCount + 1, jj + 1)
            Next jj
            ExpansionValue = Run(function_name, x_e) ' compute f(x_e)

            If ExpansionValue < ReflectionValue Then ' if f(x_e)<f(x_r) then
                For jj = 1 To ParamCount 'replace x_n+1 with x_e
                    Simplex_Matrix(ParamCount + 1, jj + 1) = x_e(jj)
                Next jj
                Simplex_Matrix(ParamCount + 1, 1) = ExpansionValue 'replace F(x_n+1) with F(x_e)
            Else
                For jj = 1 To ParamCount 'replace x_n+1 with x_r
                    Simplex_Matrix(ParamCount + 1, jj + 1) = x_r(jj)
                Next jj
                Simplex_Matrix(ParamCount + 1, 1) = ReflectionValue 'replace F(x_n+1) with F(x_r)
            End If

        ElseIf Simplex_Matrix(1, 1) < ReflectionValue And ReflectionValue < Simplex_Matrix(ParamCount, 1) Then ' else if f(x_1)<f(x_r)<f(x_n) then
            For jj = 1 To ParamCount 'replace x_n+1 with x_r
                Simplex_Matrix(ParamCount + 1, jj + 1) = x_r(jj)
            Next jj
            Simplex_Matrix(ParamCount + 1, 1) = ReflectionValue 'replace F(x_n+1) with F(x_r)

        ElseIf Simplex_Matrix(ParamCount, 1) < ReflectionValue And ReflectionValue < Simplex_Matrix(ParamCount + 1, 1) Then ' else if f(x_n)<f(x_r)<f(x_n+1) then
            For jj = 1 To ParamCount ' compute outside contraction x_0
                x_o(jj) = (1 + rho * gamma) * x_bar(jj) - rho * gamma * Simplex_Matrix(ParamCount + 1, jj + 1)
            Next jj
            OutContractionValue = Run(function_name, x_o) ' compute f(x_o)

            If OutContractionValue < ReflectionValue Then ' if f(x_o)<f(x_r) then
                For jj = 1 To ParamCount 'replace x_n+1 with x_
                    Simplex_Matrix(ParamCount + 1, jj + 1) = x_o(jj)
                Next jj
                Simplex_Matrix(ParamCount + 1, 1) = OutContractionValue 'replace F(x_n+1) with F(x_o)
            Else
                ShrinkFlag = True ' compute shrinkage later
            End If
        Else
        For jj = 1 To ParamCount ' compute inside contraction x_i
            x_i(jj) = (1 - gamma) * x_bar(jj) + gamma * Simplex_Matrix(ParamCount + 1, jj + 1)
        Next jj
        InContractionValue = Run(function_name, x_i) ' compute f(x_i)

        If InContractionValue < Simplex_Matrix(ParamCount + 1, 1) Then ' if f(x_i)<f(x_n+1) then
            For jj = 1 To ParamCount 'replace x_n+1 with x_i
                Simplex_Matrix(ParamCount + 1, jj + 1) = x_i(jj)
            Next jj
            Simplex_Matrix(ParamCount + 1, 1) = InContractionValue  'replace F(x_n+1) with F(x_i)
        Else
            ShrinkFlag = True ' compute shrinkage later
        End If
    End If

    If ShrinkFlag = True Then
        For ii = 1 To ParamCount ' compute shrinkage
            For jj = 1 To ParamCount
                Temp_x_vec(jj) = Simplex_Matrix(1, jj + 1) + sigma * (Simplex_Matrix(ii, jj + 1) - Simplex_Matrix(1, jj + 1))
                Simplex_Matrix(ii, jj + 1) = Temp_x_vec(jj)
            Next jj
            Simplex_Matrix(ii, 1) = Run(function_name, Temp_x_vec) ' compute f(x_shrink)
        Next ii
        ShrinkFlag = False
    End If

    Simplex_Matrix = BubbleSort(Simplex_Matrix)
    TotalIt = TotalIt + 1
    Loop
    
    Dim output() As Variant:    ReDim output(3, ParamCount + 1)

    output(1, 1) = "F Min Val"
    output(2, 1) = Simplex_Matrix(1, 1)
    For ii = 1 To ParamCount
        output(1, ii + 1) = "x(" & ii & ")"
        output(2, ii + 1) = Simplex_Matrix(2, ii + 1)
    Next ii
    output(3, 1) = "Iterations"
    output(3, 2) = TotalIt
    Nelder_Mead_Min = output

End Function

Function BubbleSort(ArgIn As Variant)
'x Adapted from VBA Developers GudIe with amendment to sort a matrix instead of a vector

    Dim Target_Mat As Variant: Target_Mat = ArgIn
    Dim RowCount As Integer:    RowCount = CInt(UBound(Target_Mat, 1))
    Dim ColCount As Integer:    ColCount = CInt(UBound(Target_Mat, 2))

    Dim SortedFlag As Boolean
    Dim ii As Integer
    Dim jj As Integer
    Dim bb As Integer

    Dim Temp() As Double
    ReDim Temp(ColCount) As Double

    ii = 0
    Do While (ii < RowCount) And Not SortedFlag
        SortedFlag = True
        ii = ii + 1
        For jj = 1 To RowCount - ii
            If Target_Mat(jj, 1) > Target_Mat(jj + 1, 1) Then
                For bb = 1 To ColCount
                    Temp(bb) = Target_Mat(jj + 1, bb)
                    Target_Mat(jj + 1, bb) = Target_Mat(jj, bb)
                    Target_Mat(jj, bb) = Temp(bb)
                Next bb
                SortedFlag = False
            End If
        Next jj
    Loop
    BubbleSort = Target_Mat
End Function

Let’s begin by looking at an example where we have a single feature in one dimension.  Positive class in blue cannot be linearly separated from the negative class in red.  You can clearly see that if we use a hard margin svm then our algorithm will not converge since the data is not linearly separable.  If we were to use the soft margin classifier we would be able to calibrate the model but the model performance will be terrible.

Below is another example, but this time in 2D.  Any attempts to use a linear decision boundary will fail miserably.

Feature Engineering:

For example, in the first problem presented we can change our SVM decision boundary from wx+b=0, where w and x are 1D each, to a hyperplane wz+b=0 using new feature vector z=[x,x*x].

If we plot the new feature space we can clearly see that the classes are now linearly separable.

For the second example that we have we can again increase the dimension of our feature space and see if we can linearly separate the data.  So we change our model from wx+b=0 where x=[x1,x2]  to wz+b=0 where z=[x1,x2,x1^2+x2^2].  Plotting the classes in the new 3D feature space we can see that we can indeed separate the classes with a hyperplane.

It is standard to denote feature mapping with the Greek letter Phi.  For the first example we have:

For our second example we would write:

Recall that for our soft margin classifier we derived the dual optimization as:

Subject to below two constraints:

So now that we decided to use our feature mapping function to increase the dimension of our feature space we can insert these functions into the SVM optimization function.  We get:

Next lets implement our SVM classifier just as we did in my previous posts to see if we can find a separating hyperplane in the higher dimension.

Below is my R implementation.  I first create the data and plot it.

Next we explicitly add a new feature to our data frame to increase the dimension of our feature space.

Here our Phi function would be defined as:

You can see from the 3D scatterplot that the two classes are linearly separable so let’s fit SVM to this data frame and compute our w and b coefficients.

We can see that there are 4 support vectors.

I also set this problem up in excel and solved it using Solver. The workbook link is provided later in the post.  The relevant tab is ExplicitFeatureMapping.

You can see that our results match R closely.  The difference is due to convergence criteria in R and excel.  The interesting thing I want you to note is that after solving our problem in the higher dimension we end up with the following decision boundary:

Plugging into our decision boundary equation

We get:

The key thing to note is that while we solved for a hyperplane in the higher dimension z we ended up with a non-linear decision boundary in our original feature space x.  You can see a plot of the decision boundary in the excel file.  The support vectors are highlighted in the plot.  While in this case the support vectors are close to our decision boundary they need not be when we are looking at our feature space x.  They could have been further away and not make much intuitive sense to us because they were solved for in feature space z not x.

To recap so far, we mapped our features to a higher dimension because our classes were not linearly separable.  We then searched for a linear classifier in the new feature space.  Once we find the separating hyperplane we plot it as a decision boundary in our original feature space and end up having a non-linear classifier.

Up to this point, if you have gone through all my posts on SVM, we have a method to find a large margin classifier for linear and non-linear data.  Furthermore, we have a way to allow for some violations of our margin via a tuning cost parameter.  It appears as though we have the tools to handle most practical classification problems.

Now for the bad news, there are two potential problems when we go into higher dimension feature space.  The first one is that when we increase the dimension of our features we run the risk of over fitting the model.  This is actually not a problem because we are using a model that tries to find the largest margin which controls for over fitting.  If we use cross-validation to select our tuning parameters we should end up with a model that generalizes well.  The second problem is more serious, it is that the dimensionality of the problem explodes quickly.  Suppose we have feature x1 and x2 and we wish to perform a 2-degree polynomial expansion of our features.  This means we would have the following:

This is very manageable but now think about what would happen in an image classification problem where we might be working with a 18×18 pixel image.  That gives us 324 features (one for each pixel).  If we decide that we need a non-linear SVM implementation and wanted to perform polynomial feature expansion we can see that the dimension of the feature space begins to be very large.

Kernel Trick:

There is a very cool trick to handle the problem stated above.  Notice that our data enters the algorithm only as a dot-product (highlighted in red).

What if we had a function that takes the original data and its output was dot product of the expanded features?  This will make things a lot more manageable from a computational point of view.   What I am saying is instead of explicitly mapping our features x into higher dimension z using phi and then taking the dot-product (like we did in previous example) lets use some function K(,) that takes our features x as input and returns the dot-product of z.  We can write this as:

K is called a kernel and it has to satisfy some technical conditions.  In the useful resources section I link to Mercer’s condition wikipedia entry for those who are interested.

Some of the most popular kernels are:

In the above, gamma, b, d, and sigma are tuning parameters.

Once we settle on a kernel we no longer need to explicitly map our features to higher dimension using Phi.  Instead we plug our features x into our kernel and the output is the dot-product of the higher dimensional z space.  I should mention that we don’t need to know what z space is.  In case of radial basis kernel the space has infinite dimension!

So with our chosen kernel our algorithm becomes:

The constraints remain the same.

Rewriting in matrix notation we have:

Once we find our alphas we need to recover our vector w and coefficient b so we can classify new data points.

Our classification rule for some point x star is given as:

From my previous post we derived w as a function of alphas from KKT conditions that must hold at the solution of our optimization problem.

Let’s plug this into our classification rule to get:

Now let’s remember the formula for coefficient b.  Recall from my previous post that b can be calculated from a support vector (point with a positive alpha) that has a slack variable of zero, call it point j:

Once again, plug in our definition of w:

Notice I have been highlighting the dot-products.  We again use our kernel to get the following decision boundary:

R Implementation:

I first recreated our dataset and performed a search for optimal gamma and cost parameters using a radial basis function kernel.

We can now run the SVM algorithm with C equal to 4 and gamma equal to .5.  We will later compare these results to our excel implementation.

Finally we can plot the decision boundary in the original feature space.  Please note that for some reason e1071 package in R reverses the axis in the SVM plot.

Excel Implementation:

In the workbook we have a tab with ImplicitFeatureMappig(Kernel) tab that implements this model in excel.  We first calculate kernel function values using array formulas in range C74:AE102.  Then we calculate the quadratic coefficients matrix in range C41:AE69. The objective function value is in cell D3.  The constraints are in column H.  Then we can run Solver to get the alphas.  Coefficient b is calculated in cell R5 using our second support vector.  Our results are very close to R’s output. svm_kernel_calcs_example

Hope you enjoyed.

Some Useful Resources:

• Michael Jordan’s “The Kernel Trick” in Advanced Topics in Learning & Decision Making course notes http://people.eecs.berkeley.edu/~jordan/courses/281B-spring04/lectures/lec3.pdf
• Andrew Ng notes http://cs229.stanford.edu/notes/cs229-notes3.pdf
• Lecture 15 on Kernel Methods from Yaser Abu-Mostafa at Caltech https://work.caltech.edu/telecourse.html
• Barnabas Poczos slides on Kernel Methods http://www.cs.cmu.edu/~aarti/Class/10701_Spring14/slides/kernel_methods.pdf
• R e1071 package description https://cran.r-project.org/web/packages/e1071/e1071.pdf
• Wikipedia page for Mercer’s condition that kernels need to satisfy https://en.wikipedia.org/wiki/Mercer%27s_theorem#Mercer.27s_condition.

Hard Margin Classifier Recap

We decided to use a hyperplane that separates our data perfectly with an additional requirement for that hyperplane to have the largest possible margin.  In case of two dimensions we can visualise this as a line separating class -1 from class +1.   Our formula for the separating hyperplane is given by wx+b=0.  Vector w is a weights vector that determines the slope of the hyperplane and the coefficient b determines the distance from the origin.  To derive the largest margin hyperplane we introduced supporting hyperplanes that were symmetrically placed on either side of the separating hyperplane and searched for vector w and b that maximized our margin delta.

Mathematically the problem was formulated as:

We then normalized the distance between our separating and supporting hyperplanes to have a distance of 1 on either side and made some simplifications to the formulation so that the optimization can be stated as a quadratic programming problem.

The final primal optimization problem was derived to be:

We then set up a Lagrange function and derived the dual formulation of the problem that we can pass onto the quadratic programming routine to solve for our coefficients:

Non Separable Case:

Let’s now assume that we have a dataset where the two classes cannot be linearly separated.  Most real world classification problems deal with classes that cannot be perfectly separated so this assumption is more realistic than what we dealt with earlier.  If we were to use the formulation that we derived above we would fail to find a solution because the condition y(wx+b) is equal to or greater than 1 cannot be satisfied.

To handle datasets where the classes are not linearly separable we can either use feature engineering (I will discuss this in a future post when dealing with Kernels) or we can relax the constraint that is causing us grief.

To relax the constraint that data is linearly separable we can introduce a new slack variable.  The constraint will now look like the following:

To visualize this I have added two new data points to our graph and circled them.  You can see that if we use the blue hyperplane as our classifier we would be making two errors.  The slack variables measure the extent to which the margin, y(wx+b)=1 is violated.

We can have 3 cases for any particular data point.

• Correct classification and data point does not violate the margin. In this case slack variable is equal to zero
• Correct classification but data point violates the margin. In this case slack variable is between zero and one.  This will happen when a data point is on the correct side of the separating hyperplane but is inside the margin that is defined as wx+b=1.
• Incorrect classification. This is the case that is depicted in the graph above.  The data point is on the wrong side of the separating hyperplane and y(wx+b)<0 and therefore the slack variable is greater than one.

Now that decided to introduce the slack variables we need to modify our objective function.  To balance any violations of our margin we want to penalize our objective function.  To do so we will introduce a total cost measurement.  This will capture the total number of violations we make.  I should point out that this total cost doesn’t measure the classification error rate; we are only measuring the violation of our margin.  A point can be on the wrong side of the margin but still classified correctly.  Let’s define total cost as:

We scale the sum of our slack variables by a tuning parameter C.  The value of C is not solved for in our optimization routine but is selected via cross validation.

So remembering that for a linearly separable dataset our optimization problem was:

And now with the introduction of the slack variable that allows us to fit a separating hyperplane to a non separable dataset we have the following primal optimization problem:

So we are trading off total cost versus width of the margin.  Tuning parameter C controls this trade-off.  If we were to select some large value for C then we will penalize the model for any violations of the margin. In the limit as C is increased we are left with the case of a hard margin classifier.

To solve this problem we need to set up a Lagrange function which captures our amended objective function and the two constraints:

We are minimizing our objective function with respect to w, b, and the slack variable but maximizing with respect to alpha and beta.

The KKT conditions for the solution are given as:

Prior to differentiating the Lagrange function lets group variables together to make things more clear

Now differentiating with respect to the vector w we only use the blue terms to get:

By setting the derivative to zero we have:

Now we differentiate with respect to coefficient b (the green term) to get and set the results equal to zero:

And finally we can take the derivative of the Lagrange function with respect to the slack variable:

And as with the previous derivatives we set it to zero:

We will be using the KKT conditions and formulating the dual optimization problem to achieve some remarkable results so let’s rewrite our Lagrange as:

Notice that the green and red terms equal to zero at the solution due to our KKT conditions.

So we have:

We can now substitute for the value of w what we derived from one of the KKT condition above:

And here is the very cool result; we end up with an identical set up as in the hard margin case from Part 1

The dual problem is given by:

Subject to the following constraints:

The last constraint is due to the fact that at the solution KKT condition C-alpha-beta =0 must hold.  Since alpha and beta must be equal to or greater than zero then C must be greater than or equal to alpha.

So notice what we were able to show.  We introduced slack variables to allow us to handle non-separable data but after deriving the dual formulation of the problem we are left with something that is very similar to the hard margin problem with the exception of an additional constraint that alphas must be smaller than some tuning parameter C that a user selects.

Once we use a solver to calculate our alphas we need to work out what our w vector and parameter b is to use the classifier.  Again, using KKT condition we were able to derive the value of w given alphas as:

Once we calculate vector w we can use the fact that for our support vectors (data points that have non-zero alphas) and that are not violating the margin (slack variable is equal to zero) we have:

We solve for b to get:

Once we have w and b we can use the classifier by simply plugging in some observation x and using the rule:

Excel Example:

In the workbook I set up the dual formulation of the quadratic programming optimization.  We are trying to maximize cell D4 by changing values of alpha.  Once we solve for aphas we can calculate w and b in cell M2:M4.  With our coefficients in hand we can assign a class label to a new data point in cells P4 and Q4.  We can control the cost parameter C in cell D3. svm_calcs_example

Sanity Check:

As in Part 1 we can run a quick check to make sure our implementation matches R.

Our results match so we can have some confidence that our spreadsheet is set up correctly.

Hope you enjoyed.

Some Useful Resources:

• Idiot’s Guide to SVM from MIT http://web.mit.edu/6.034/wwwbob/svm-notes-long-08.pdf
• Andrew Ng notes http://cs229.stanford.edu/notes/cs229-notes3.pdf
• Lecture 15 on Kernel Methods from Yaser Abu-Mostafa at Caltech https://work.caltech.edu/telecourse.html
• Christopher Burges “A Tutorial on Support Vector Machines for Pattern Recognition” http://research.microsoft.com/pubs/67119/svmtutorial.pdf
• Practical Mathematical Optimization by Jan A Snyman is excellent http://www.springer.com/us/book/9780387243481

I decided to check out some other online content on SVM and have a stab at presenting the material in a style that is friendly to newbies.  This will be a three part series with today’s post covering the linearly separable hard margin formulation of SVM.  In a second post I will discuss a soft margin SVM. In a final post I will discuss what is commonly referred to as the kernel trick to handle non linear decision boundaries.  I borrow heavily from the authors that I will mention at the end of the post so to add some original content I implement SVM in an excel spreadsheet using Solver.  I find that excel is a great pedagogical tool and hope you agree with me and find the toy model interesting.

The Goal:

In this post I will work with an example in two dimensions so that we can easily visualize the geometry of our classifier.  All the points discussed apply equally in higher dimensions.

We will work with a binary classification problem.  Our features space will consist of 12 observations of x1 and x2 pairs.  In general we will have N observations and P will denote the number of features.  Each data point is labeled as a +1 or a -1.  In the plot I have + for +1 and – for -1.

What a linear classifier attempts to accomplish is to split the feature space into two half spaces by placing a hyperplane between the data points.  This hyperplane will be our decision boundary.  All points on one side of the plane will belong to class C1 and all points on the other side of the plane will belong to the second class C2.

Unfortunately there are many ways in which we can place a hyperplane to divide the data.  Below is an example of two candidate hyperplanes for our data sample.

Both H1 and H2 separate the data without making any classification errors.  Points above the line are classified as +1 and points below the line are classified as -1.   A solution to this problem of finding an ‘optimal’ separating hyperplane is to introduce the concept of a margin.  We first select a hyperplane (decision boundary) and then extend two parallel hyperplanes on either side of our decision boundary.  We widen the parallel hyperplanes until we touch one or more of our data points on either side but without allowing any of the points to penetrate the boundaries.  The distance between these parallel hyperplane boundaries is called a margin.  The decision boundary is equidistant from these parallel hyperplanes so that we do not bias our decision in favour of any particular class.  In below plot I show these margins for each decision boundary.  For H1 the parallel hyperplanes are plotted in green.  For H2 the parallel hyperplanes are plotted in orange.  The margins are plotted as dashed lines.  You can see that the margin is larger for H2.  You can think of the margin as the width of a road.  We try to build a straight road that is as wide as possible between the data.  The separating hyperplane is the dividing lane on this road.  We are trying to find the placement of the dividing lane in such a way as to build the widest road that just touches any of our data points.

Defining a hyperplane:

So now that we have our goal set out for us we need to formulate it as a mathematical model.  As we plow along I will take lots of detours to introduce the algebra and geometry.  With each detour I will try to keep us focused on why the concept is important for our problem.

With that in mind let’s begin with a definition of a vector.  A vector can be thought of as a directed line segment.  A vector has a displacement from one point P to another point Q.  A starting point of a vector is called a tail and the ending point is called a head.  If a vector starts from the origin it is called a position vector.  Below is a plot of a two position vectors in two dimensions.

Vectors do not need to have a tail at the origin.  Two vectors are equal if they have the same length and the same direction.  Below we have vector v from O to A and vector w from B to C.

Notice that vector v has displacement of [3,2]. Vector w displacement is also [3,2] and is given simply as the difference of its two components [3,2] = [8-5,3-1] or C-B.

We will use vectors throughout this post.  They are critical to our definition of a hyperplane and we will also use operations on vectors to measure the margin.  Before we move on I want to mention the dot product.  For vectors u and v we define a dot product as:

With a definition of a dot product we can now specify an equation for our hyperplane:

In our two dimensional example, x is a vector of [x1,x2].  Vector w and scalar b are the coefficients of our hyperplane.  If for example w=[1,2] and b=3 then our hyperplane is given as:

We can rewrite this in the more familiar y=mx+b form by rearanging some terms.

Notice that the slope of the line is given by w1/w2 and the intercept term is given by b/w2.

Let’s plot this hyperplane and also include the vector w:

The important thing to note here is that the vector w is perpendicular to our hyperplane.  This means that vector w is at a 90 degree angle to our hyperplane.  This is a very important property, referred to as orthogonality, which we will exploit later when calculating our margin.  The vector w is called the normal vector.

Another important property that will be relevant in our discussion is that multiplying vector w and b by a constant (scalar multiplication) does not change our hyperplane but does extend our normal vector in the original direction.  To see this lets multiply w and b by 2 and plot.  We get:

So now we know that our separating plane will be defined as:

And our decision rule for the model will be:

That is, we will assign observation xi to class +1 if wx+b is >0 and assign xi to class -1 if wx+b <0.

Defining the Margin:

Lets again look at a plot of what we are trying to measure, and then I will introduce the bare minimum mathematical definitions for us to be able to compute what we need.  So here goes:

What we need is the distance from a point to our hyperplane.  In the above graph we are interested in the distance delta that is marked as a dashed orange line.  To do so we need to first define the length or magnitude of a vector.   Euclidian length (or norm) of a vector v is defined as the square root of the dot product of v:

Above result is given by the Pythagoras’ Theorem.  In two dimensions we can see this more clearly. Below I plot vector v=[a,b] and calculate its length from the origin.

Once we have the vector length in hand we can normalize a vector to have unit length.  We simply divide each entry in the vector by the norm.  For example assume we have a vector v= [4,3]. Then unit vector u is:

And finally, let’s define a measure of distance between two vectors.  We can use our measure of length on the difference between two vectors to get its Euclidian distance.

For example, if we have vectors u and v, the distance between them is calculated as:

We now have all the tools we need to calculate the delta value.

In below plot we need to derive the value of delta (dashed orange line).

This value of delta is ½ the margin for our decision boundary.  In Euclidean geometry, the shortest distance from a point to a hyperplane is one that is orthogonal to the hyperplane (perpendicular to the line in 2 dimensions). In the plot, delta equals the difference between vector A and vector B.  We know that delta is in the opposite direction of our normal vector (marked in red).  We can use this fact and first compute the pure direction of our normal vector by dividing it by its norm and then multiply it by delta to gives us a vector of length delta and pointing in the direction of the norm.  We then subtract this value from A to arrive at point B on our hyperplane.

Let’s pause and understand what we did.  We are trying to find a point on the plane that is closest to one of our observations xi.  In the plot we have labeled the point we are trying to compute as B and the data point from where delta is measured is labeled A.  Remember that A is one of our actual observations xi.   We then multiply the length delta by the direction of our normal vector and then subtracted this quantity from our observation xi to end up on the hyperplane at point B.  Once we have point B we can use the fact that any point on our hyperplane produces a value of zero.  Therefore, we can plug in B into the formula for our hyperplane and solve for delta.

Remember our definition of norm:

Plugging that in and solving for delta we have:

Notice that in the calculations so far we worked with xi that is on the positive side of the decision boundary.  We need to take the absolute value of our delta calculations because when computing the delta from the negative side of the decision boundary we will have a negative value.  However, instead of taking the absolute value we can simply multiply our delta by the class label.  So our final formula becomes:

When we are below the decision boundary the value in the bracket is negative and it is multiplied by -1 (our class label).  So we are now dealing with strictly positive deltas.  The above formula is commonly referred to as the geometric margin.

So we now have a formula to compute the distance between any of our observations to our decision boundary.  Lets now remember that delta is defined here as just a single lane of the road (ie, it is half the margin).  Therefore, 2*delta is our margin for any point.  The definition of a margin for the entire sample is simply the minimum margin of all the possible margins in the sample.  And it is important to remember that we are trying to maximize the sample margin.

Therefore, our mathematical formulation for the hard margin classifier that we defined at the beginning is:

Optimization Set Up:

Now that we have a mathematical formulation for a maximum margin classifier we have to start thinking about how do we find the w and b that achieve the maximum.  The rest of the post is about simplifying this maximization problem into something that we know we can solve.  To begin with let’s recognize that our definition for the delta as the minimum delta in our sample can be restated as:

Our objective function remains unchanged and we only changed the way we represent the delta of a sample.  This equation reads as, maximum of margin with respect to w and b where for all xi and yi in our sample the geometric margin is at least equal to delta.

Also, we need to come back to the important point I made earlier about length of w not having any impact on our separating hyperplane.  Remember that if we rescale our normal vector the slope of the hypeplane does not change.  So for convenience let’s set norm of w equal to:

I realize that this is not an obvious choice but if we now solve for delta to get one over length of w and substitute into our maximization formula we can see that this transformed our objective function and the constraints:

On the right side in the constraints equations please notice we can cancel out the norm of w to leave us with linear constraints.  Linear constraints are much easier to work with so we should be happy with our accomplishments so far.

A few final changes are needed to change this maximization problem into something we can easily tackle.  The major problem now is the objective function itself.  The one over w norm is not a nice equation that is easily maximized.  What we can choose to do instead is to minimize the reciprocal of the objective function.

Finally, we can get rid of the square root in the objective function so that we can work with a very well known quadratic programming problem.  We can do this because the square of a function is a monotonic transformation.   In a quadratic programming problem we work with a quadratic objective function subject to linear constraints.  There are many efficient algorithms in most software packages that can solve these problems.

The above problem is what is known as the Primal Form of the optimization problem.

Excel Example 1:

We can now set up a toy example in excel. Below is an example for our data set.

I have included a screenshot with the formulas.  The objective function is in cell C2.  We can add constraints in Solver and then hit <Solve> button.  As can be seen we have 3 data points for which the constraint in the minimization is binding (ie geometric margin is equal to one).  In the plot I highlight these points.  These are called the supporting vectors for our margin classifier.

The Dual Set Up:

To solve the constrained optimization that we have derived so far we can set up a Lagrange function and derive the solution analytically.  Let’s restate the minimization we are trying to compute:

We can rewrite the constraint as:

We then multiply all of the constraints by Lagrange multipliers (alphas) and subtract from our objective function:

Notice that we are minimizing with respect to w and b but are maximizing with respect to alphas.

As I mentioned, the quadratic optimization problem is a well studied problem.  A quadratic objective function with linear inequality constraints must satisfy some conditions at the solution.  These conditions are known as KKT conditions and are named after Karush, Kuhn and Tucker who independently derived them.

The KKT conditions are:

First lets expand the Lagrangian to isolate the w and b terms:

Now we can take the derivative of our Lagrange function with respect to the vector w to get:

And by setting the result equal to zero and solving for w we get:

Now we can move on to the derivative of the Lagrange formula with respect to b.  We have:

And we again set the result equal to zero:

Now let’s substitute these results into the Lagrange function:

It is important to note that we are now trying to maximize this Lagrange function since it is only a function of alphas.  This is the so-called Dual formulation of the optimization problem.

This can be rewritten in matrix notation as:

Once we solve for alphas we need to be able to use our classifier.  That means we need the coefficients of our hyperplane.  To get back our coefficients w we can use below equation that we derived earlier:

Once we have our w we can use the fact that for any of the support vectors we have:

And we can solve this for our b coefficient:

Now that we have our hyperplane coefficients we can use our classifier to assign any observation xi by plugging it into below equation:

Excel Example 2:

We can set up the maximization problem in a spreadsheet and use Solver to calculate alphas.  Below are the results.  You can see that we have only three alphas that are non zero and they match the support vectors we derived in the Primal example.

I have attached the sheet svm_calcs_example

Sanity Check:

We can perform a quick spot check in R to make sure our results are correct.

Fortunately we got results that are the same as our excel implementation.

Hope you enjoyed the post. Below are some useful resources.

Useful Resources:

• Idiot’s Guide to SVM from MIT http://web.mit.edu/6.034/wwwbob/svm-notes-long-08.pdf
• Andrew Ng notes http://cs229.stanford.edu/notes/cs229-notes3.pdf
• Lecture 14 on SVM from Yaser Abu-Mostafa at Caltech http://work.caltech.edu/lectures.html
• Thorsten Joachims SVM lecture from Cornell http://svmlight.joachims.org/
• Practical Mathematical Optimization by Jan A Snyman is excellent http://www.springer.com/us/book/9780387243481
• Data Mining by Nong Ye is a good book that has a worked example for the AND function https://www.crcpress.com/Data-Mining-Theories-Algorithms-and-Examples/Ye/p/book/9781439808382
• I used the wife’s Linear Algebra and Its Applications for Algebra review but any textbook will do. https://www.amazon.com/Linear-Algebra-Applications-Updated-CD-ROM/dp/0321287134/ref=sr_1_7?s=books&ie=UTF8&qid=1465859903&sr=1-7

Dataset Description

• Private Nonfarm Payrolls. The most important statistic in the market and generates the most volatility.  Released on first Friday of a month and measures the monthly change in the number of employees on business payrolls. This is the target variable in the models.  I also include lagged values to capture autocorrelation.
• Temporary Help Services. This statistic is released on the same day as the Nonfarm payrolls data and captures shifts in hiring of temporary workers.  Increase in temp hiring is potentially a leading indicator for the labour market and is lagged in my models. I use up to 3 historical monthly readings.
• ADP Private Nonfarm. Automatic Data Processing Inc is a leading payroll related services provider. Data from 3650 of their private clients is used to construct a forecast of the official NFP release.  The model used to generate the NFP forecast is built in collaboration with Moody’s Analytics.  This is a timely report since it precedes NFP by two days.
• Initial Claims 4 week average. A moving average of the weekly initial unemployment claims.
• Continuing Claims MoM change. Monthly change in the number of people that have filed jobless claims with the appropriate government labour office.
• Manufacturing PMI Employment Index. A diffusion index of a monthly survey by ISM of supply management professionals.  Readings above 50 means factories are hiring new workers.
• Non-Manufacturing PMI Employment Index. Similar to the manufacturing index but surveys the service industry.
• Conference Board Labour Differential Index. This is a report that surveys approximately 3000 households about their confidence about the labour market.

Model List:

Below is a sample list of some popular regression modeling approaches in machine learning which I will consider.  R’s caret package has a function train that can handle all these models and much more.  More on this package at http://topepo.github.io/caret/

• Ridge Regression [100 models].   
  • Tuning parameter: Labda
• Support Vector Regression with Linear Kernel [100 models].
  • Tuning parameter: Cost
• Support Vector Regression with Polynomial Kernel [400 models].
  • Tuning parameter: Polynomial Degree, Cost, Scale
• Support Vector Regression with Radial Basis Function Kernel [500 models].
  • Tuning parameter: Cost, Sigma
• Stochastic Gradient Boosted Trees [684 models]
  • Tuning parameter: Number of trees, Maximum depth of tree, Shrinkage
• Cubist [36 models]
  • Tuning parameter: Number of trees, Maximum depth of tree, Shrinkage
• Neural Network [50 models]
  • Tuning parameter: Number of nodes in a hidden layer, Weight decay
• Nearest Neighbor Method [30 models]
  • Tuning parameter: number of clusters

In total I will consider 1900 models. For each class of models repeated 10-fold cross validation will be performed with 5 repeats to select optimal tuning parameters.  Fifty monthly observations will be left out for testing model performance.  In total I have 131 monthly observations for the period starting Jun 2005 and ending Apr 2016.

my.data<-read.csv('nfp2.csv')
(n<-nrow(my.data))

train.data<-my.data[51:nrow(my.data),] # data for model tuning
test.data<-my.data[1:50,] # data for model generalization

require(caret) # caret package for train function
train.cntrl<-trainControl(method = 'repeatedcv', number = 10, repeats = 5)

In today’s blog I will present the R code for Ridge Regression, SVM with linear and polynomial kernel.  Caret’s github page has the details about all other models and readers should be able to easily replicate my results.

Ridge Regression Results:

We can estimate the model using below code

####################################
# Ridge Regression Model
####################################
set.seed(81)

ridge.grid<- data.frame(.lambda = 10^seq(1, -5, length = 100))

ridge.fit<-train(train.data$NonFarm~.,
data=train.data,
method='ridge',
tuneGrid=ridge.grid,
trControl=train.cntrl,
preProc =c('center','scale'))

plot(ridge.fit) # plot CV profile
ridge.fit$results[ridge.fit$bestTune[[1]]==ridge.fit$results[,1],] # return RMSE for optimal lambda

ridge.test.predict<-predict(ridge.fit, newdata=test.data) # test data fit
(ridge.test.rmse<-sqrt(sum((test.data$NonFarm - ridge.test.predict)^2)/nrow(test.data))) # test RMSE

The cross validation RMSE vs lambda parameter profile is presented below:

The optimal lambda parameter of .0142 results in RMSE of 75.56 and Rsquared of 84%.

SVM- Linear Kernel Results:

To estimate Support Vector Regression model with a linear kernel we can use below code:

####################################
# svm linear kernel
####################################

set.seed(81)

svm.linear.grid<- data.frame(.C = 1:100)

svm.linear.fit<- train(train.data$NonFarm~.,
data=train.data,
method = 'svmLinear',
trControl = train.cntrl,
tuneGrid=svm.linear.grid,
preProc =c('center','scale'))

The cross validation plot is below:

The optimal cost parameter of 1 results in RMSE of 74.95 and Rsquared of 85%.

SVM- Polynomial Kernel Results:

To estimate SVM with polynomial kernel we set up a grid of tuning parameters and use the train function as in previous examples.

####################################
# svm polynomial kernel
####################################

set.seed(81)

svm.poly.grid<- expand.grid(.degree =c(2,3),
.scale=10^seq(2, -5, length = 20),
.C = 2^seq(-2, 7, 1))

CV RMSE profile:

The optimal model has degree of 3, scale of .0002976, and a C value of 128. The reported RMSE is 74.74 and R squared of 85%

SVM- Radial Basis Function Kernel Results:

With radial basis kernel we have two tuning parameters

####################################
# svm radial basis
####################################

set.seed(81)

svm.rbasis.grid<- expand.grid(.C = 2^seq(-2, 7, 1),
.sigma=seq(1 , 50, by=1 )/100)

The optimal cost parameter is 8 and vol of .01.  Best RMSE is 75.13 and Rsquared is 85%

Boosted Trees Results:

Gradient Boosting Machines are all the rage so lets see how they fit our data.  We can use the gbm package in R. Lets set up a tuning grid and check cross validation RMSE.

####################################
# gbm
####################################

set.seed(81)

gbm.grid<- expand.grid(.interaction.depth = seq(1, 7, by = 2),
.n.trees = seq(100, 1000, by = 50),
.shrinkage = c(0.01, .05, 0.1),
.n.minobsinnode=c(5,10,20))

Optimal tuning values for the model were 550 trees with a depth of 1, a shrinkage parameter of .01. and minimum number of observations in a node of 5.  RMSE is 83 and Rsquared if 82%.

Cubist Results:

Optimal tuning values of 0 neighbors and 5 commitees returned RMSE of 76.98 and Rsquared of 84%.

####################################
# cubist
####################################
set.seed(81)

cubist.grid<- expand.grid(.committees = c(1, 5, 10, 50, 75, 100),
.neighbors = c(0, 1, 3, 5, 7, 9))

Neural Network Results:

A neural network with a single hidden unit and weight decay of .25 resulted in RMSE of 91 and Rsquared of 79%

####################################
# nnet
####################################
set.seed(81)

nnet.grid<- expand.grid(.decay = c(0.01, .05, .1, .25, .5),
.size = 1:10)

Nearest Neighbor Results:

This is the last model that I fitted.  The optimal number of clusters was 6 which resulted in a cross validation RMSE estimate of 88.9 and Rsquared of 78%.

####################################
# knn
####################################
set.seed(81)

knn.grid<- data.frame(.k = 1:30)

Model Comparison:

We can also perform inter model comparison in R.  Because we used the same repeated cross validation specification for all the models and we set the same seed prior to training we have paired samples ready for comparison.  With 10-fold cross validation and 5 repeats we have 50 resamples for each model.  We can use caret’s resamples and diff.resamples functions to check if the difference in model  performance metrics are statistically significant.

#######################################
# Check model diff
#######################################
resmpl <- resamples(list(Ridge = ridge.fit,
SVM_lin = svm.linear.fit,
SVM_poly=svm.poly.fit,
SVM_rbf=svm.rbasis.fit,
GBM=gbm.fit,
Cubist=cubist.fit,
NeuralNet=nnet.fit,
KNN=knn.fit))

model.diff<-diff(resmpl)
summary(model.diff)

I exported the data into excel and formatted the table for easy comparison.

The results suggest that GBM, Neural Network and k-nearest neighbor methods are inferior.

Finally, we can check how well the models generalize by checking RMSE on our test data.  Below is a table with the first 6 columns reporting summary statistics of the RMSE metric that we obtained on the cross validation run.  The last column reports the RMSE on our test data. Unsurprisingly the cross validation RMSE is more optimistic than the test data RMSE.  Since I used data from Apr 2001 to May 2005 as the test data there is some hope that this is the upper bound on the true RMSE.  This is possible since the test data is older and we used more recent (perhaps more relevant) data in the model calibration.

Models vs Market:

This was a simple exercise in model selection and it produced modest results.  We live in a data rich world and I used a minuscule number of features compared to the data that is available.  I omitted regional labour market indexes like Empire State Manufacturing Index, PhillyFed Index, all of the detailed data that is released by BLS as part of its employment report.  Non conventional sources of data can also be included for consideration.  Google trend data on key search terms is one example.

Having said all that it is time to check how the models compare to the market consensus estimate for Nonfarm Payrolls.  Since Jun 2005 to present the RMSE of the bloomberg consensus estimate has been 62k!  This reveals that the market forecast is based on a much more rich information set than the one considered here.

Useful Resources:

• I have mentioned this in many posts, Applied Predictive Modeling is a MUST. All of the models I mentioned in this post are in that book and MUCH MUCH MUCH more http://appliedpredictivemodeling.com/
• Caret github page https://topepo.github.io/caret/
• Slides on Cubist model https://cran.r-project.org/web/packages/Cubist/vignettes/cubist.pdf
• Good book on economic data releases https://www.amazon.ca/Secrets-Economic-Indicators-Investment-Opportunities/dp/0132932075

I will use Fed’s Apr 2016 press release webpage as an example.  Our goal is to import the content of the page into R and then to parse the data so that we are left with just the text organized in some efficient data structure.

We first need to load the RCurl package and use getURL function to request and retrieve the content from a web server.  The returned object is a vector of length one that contains the retrieved information.

Our aim is to extract the main text of the page and discard all other data that was returned to our request.  To achieve this we first create a tree structure from the Apr2016.request object using htmlTreeParse function from the XML package.

If you inspect the Apr2016.tree content you will notice that the main text of the webpage is contained in the paragraph <p> tag.

We can use xpathApply function to extract the content of each paragraph.  To return a vector instead of a list we will unlist the output.

We now have the content of all the paragraphs of the Fed statement in a single R object.  To perform text mining we can use R’s tm package.

We first need to create a text corpus which is just a structured set of documents.  For that purpose we can use tm’s Corpus function.  To create a corpus object from a vector of text we should use VectorSource function.  VectorSource function interprets each element of a text vector as an individual document.  So far, our Apr2016.tree.parse vector has a paragraph as its element so we need to combine each paragraph into a vector of length 1. I use a loop and omit the first paragraph (the date of the statement) and the last paragraph (link to Fed’s policy implementation notes).

Finally we should remove \r (carriage return) and \n (new line) strings from the text and then create a corpus.

Now we can use tm’s text mining functions to change all words to lower cases, remove punctuation, strip out white space, and to remove stop words.  Stop words are commonly used words that are often removed from a corpus before any analysis is performed.  Typical stop words are ‘me’, ‘what’, ’were’, ‘who’, ‘is’, etc.

Once we have a corpus that is prepared for analysis we can change the data structure to a term document matrix.  A term document matrix has a row for each unique word and the value in each row is the number of times that particular word appears in a document.  Each column belongs to a particular document in a corpus.  In our case we have a single document so we have one column.

And that is all there is to it.  We now have the Fed statement document ready for text mining.  To conclude the post I will simply create a wordcloud of the statement just because I feel like today’s entry lacks colour. The example is taken from Nigel Lewis’ book on data visualization.

Useful Resource:

• Tutorial on tm package https://cran.r-project.org/web/packages/tm/vignettes/tm.pdf
• I use http://stackoverflow.com for all my coding needs. Great community. If you have a question about R, most likely its been answered on this forum.
• LSE “Text mining for central banks” is an interesting piece on text mining http://eprints.lse.ac.uk/62548/1/Schonhardt-Bailey_text%20mining%20handbook.pdf
• Excellent book on data visualization where I took the wordcloud example http://www.amazon.ca/Visualizing-Complex-Data-Using-R-ebook/dp/B00I7I0VPQ?ie=UTF8&*Version*=1&*entries*=0

I have previously discussed penalized regression models in the context of ridge regression and the LASSO model.  These two models are special cases of the elastic net model.  Recall that in Ridge regression we included an L2 penalty term in our sum of squared errors loss function which we attempt to minimize to estimate our theta parameters:

The lambda term is a hyper parameter and is estimated using cross validation.

A higher value of lambda shrinks the model’s coefficients towards zero.  The particular choice of using an L2 penalty term means that our estimated theta coefficients approach zero but do not equal zero exactly.

In my post on Ridge regression I presented a plot of coefficient paths as a function of lambda.

Penalized linear regression models try to balance the bias variance trade-off by imposing increased bias in exchange for a (hopefully) larger decrease in the variance of the model.

The LASSO model uses an L1 penalty term in the loss function we are trying to minimize:

The lambda parameter serves the same purpose as in Ridge regression but with an added property that some of the theta parameters will be set exactly to zero.

Elastic Net:

The elastic net model combines the L1 and L2 penalty terms:

Here we have a parameter alpha that blends the two penalty terms together.  When alpha equals 0 we get Ridge regression.  If alpha is set to 1 then we have the LASSO model.  The lambda parameter is the shrinkage coefficient.

R implementation:

To estimate the model in R we can use the glmnet package that has elastic net model implementation.  In glmnet we can perform cross validation to find the lambda parameter that returns the smallest possible root mean squared error statistic for a selected alpha parameter.  This approach is useful when we decide apriori on what alpha we want to use.  If we have resolved to use Ridge regression we can perform cross validation to find optimal lambda while keeping alpha set to 0.  Alternatively, if we wish to find the optimal lambda for the LASSO model we would set the alpha parameter equal to 1.

In our case we want to find the optimal lambda and alpha jointly. For that we will need to use the caret package.  Using the train function in the caret package we can set up a grid of alpha and lambda values and perform cross validation to find the optimal parameter values.

Below is an example using Hitters dataset from ISLR package.

Practical Model:

If I haven’t lost you with the baseball dataset I want to show you a more practical model.  Many brokers produce toy regression models as a sanity check of valuations in FX.  Below is one example produced by a shop that is particularly fond of these models.

We can use elastic net to estimate similar models that will generalize better (we hope).  In addition to using elastic net, I like to use interaction variables.  Let me digress for a minute to explain interaction variables in the context of an FX fair value model.  Usually an FX pair is regressed on rate spreads, perhaps relative equity index valuation, and let’s say the VIX index to capture risk aversion.  A model would be estimated as:

We can include interaction terms as follows:

I realize that this seems odd but we can group the variables together and then factor to get:

Notice that what we end up with are theta coefficients that are dependent on the features.  For example, the impact of rates spread is dependent on the level of relative equity valuation and the VIX.  By introducing interaction terms we can pick up interesting dynamics that usual linear regression models cannot.  If such theta-feature dependencies exist, and are not swallowed by noise, we should be able to pick them us using our elastic net model.

Elastic Net FX Model:

With a little bit of effort we can shoehorn elastic net models into excel and have live models for select FX pairs.  Below is an example of a model for EURAUD.  I cannot stress enough the importance of including out of sample test data to check how well our model generalizes.  If the model does a decent job of fitting the test (out of sample) data then we can place more confidence in our model.

There are instances where our model does well on training data but does not generalize well.  For example, AUDJPY model fits well in-sample but does a poor job out of sample.

A Warning:

I want to end on a cautionary note.  These models are useless for systematic trading.  They backtest poorly.  I believe they are extremely useful for scenario analysis or as input into discretionary trade selection but should not be used to trade FX just because the market deviates from the model.  This should be obvious to everyone.  If we make a generous assumption and assume our model captures a real relationship between the target variable and the features, we still cannot trade that signal using the FX pair.  If we observe a deviation of the market price from our model price there is nothing that forces the FX market to clear that difference.  It is just as plausible for the rate spread to realign the model price to the market value of the FX pair.

Some Useful Resources:

-I keep telling you to pick up Applied Predictive Modeling by Max Kuhn and Kjell Johnson. Its well worth the money.  http://appliedpredictivemodeling.com/

-An introduction to Statistical Learning is another great book.  Some of the authors built the glmnet package. http://www-bcf.usc.edu/~gareth/ISL/

-Glmnet vignette https://web.stanford.edu/~hastie/glmnet/glmnet_alpha.html#lin

– My post on Ridge regression https://quantmacro.wordpress.com/2015/12/11/ridge-regression-in-excelvba/

– My post on LASSO model https://quantmacro.wordpress.com/2016/01/03/lasso-regression-in-vba/