Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
I have signed up for Attilio Meucci’s ARPM Bootcamp next month (July 13-18) in NYC http://www.symmys.com/arpm-bootcamp, and need to do quite a bit of prep as it’s going to be a deep-dive…
The Advanced Risk and Portfolio Management Bootcamp provides in-depth understanding of buy-side modeling from the foundations to the most advanced statistical and optimization techniques, in 6 intensive days of theory and MATLAB live examples and exercises:
- Market modeling: random walk, ARMA, GARCH, Levy, long memory, stochastic volatility
- Multivariate statistics: non-parametric, non-normal MLE, shrinkage, robust, Bayesian estimation; copula/marginal factorization; location-dispersion ellipsoid
- Factor modeling: theory and pitfalls of time-series and cross-sectional factor models, CAPM, APT, principal components analysis, random matrix theory
- Pricing: full evaluation, Greeks, stress-matrix interpolation; analytical, Monte Carlo, historical
- Risk analysis: diversification, stochastic dominance, expected utility, Sharpe ratio, Omega, Kappa, Sortino, value at risk, expected shortfall, coherent and spectral measures
- Portfolio construction: robust/SOCP optimization, shrinkage/Bayesian allocations, Black-Litterman and beyond; transaction costs, liquidity, market impact; statistical arbitrage; convex/concave dynamic strategies, CPPI, delta-replication
So I thought I would delve into one of his many interesting papers, “Neither ”Normal" not “Lognormal”: Modeling Interest Rates Across all Regimes“, co-authored by Angela Loregian. The paper and Matlab code can be downloaded at http://www.symmys.com/node/601.
However, since I’m still an R fan I thought it would be interesting to port his code to R (although there is an R package, “Meucci”, under development, it failed to build when I looked). Here is my port of the code and the resulting chart showing JGB rates, log-rates, and shadow rates (derived from the inverse call transformation):
################################################################################
# load packaages #
################################################################################
library(RCurl)
library(xts)
library(matlab)
library(pracma)
library(reshape2)
library(ggplot2)
library(gridExtra)
################################################################################
# get JGB data from file #
################################################################################
csvfile = getURLContent(
"https://docs.google.com/uc?export=download&id=0B4oNodML7SgSZGdVQXYxVWRwMW8",
ssl.verifypeer = FALSE, followlocation = TRUE, binary = FALSE)
JGB <- read.zoo(textConnection(csvfile), format = "%Y-%m-%d", index = 1,
header = TRUE, tz = "UTC", sep = ",")
JGB <- as.xts(JGB)
# process data
tau <- as.numeric(gsub("Y", "", names(JGB)[c 1="2," 2="3," 3="5," 4="7," 5="10," 6="11," 7="12)" language="(1,"][/c]))
# rates
y <- coredata(JGB[, c(1, 2, 3, 5, 7, 10, 11, 12)])
Date <- index(JGB)
TimeStep <- 5 # daily (1) / weekly (6) observations
date <- Date[seq(1, NROW(Date), TimeStep)]
y <- t(y[seq(1, NROW(Date), TimeStep), ])
################################################################################
# Inverse Call Transformation function #
################################################################################
# this function computes the Inverse Call Transformation and returns shadow
# rates see A. Meucci, A. Loregian - "Neither "Normal" not "Lognormal": Modeling
# Interest Rates Across all Regimes" to appear (2013)
# Last version of code and article available at http://symmys.com/node/601
#
# INPUT
# rates: matrix containing the time series of the rates to be transformed;
# size(rates)= lenght(tau) x t_, where t_ is the length of the time series
# tau: vector containing the times to maturity corresponding to the rows of the
# rates matrix
# eta, zeta: Inverse-call transformation parameters (the smoothing parameter s
# is obtained as s=eta*exp(zeta*tau))
#
# OUTPUT
# x: shadow rates, computed from rates via inverse-call transformation
################################################################################
InverseCallTransformation <- function(rates, tau, eta, zeta) {
t_=size(rates,2);
x=zeros(size(rates)[1], size(rates)[2]);
s=eta*exp(zeta*tau);
# Bachelier call pricing function
BachelierCallPrice <- function(x,s) {
# Call function (zero-strike call option price profile according to the
# Bachelier pricing function)
# s: smoothing parameter
C = x * pnorm(x/s) + s * dnorm(x/s);
}
call_fit <- function(tmpX,y,sigma) {
c=BachelierCallPrice(tmpX,sigma);
F= y - c;
}
for (v in 1:length(tau)) {
# inverse call transformation
# Optimization options.
opts <- list(tau = 0.01, tolx = 1e-10, tolg = 1e-10, maxeval = 7*500)
x0=0; #initialization
for (t in 1:t_) {
# fitting inverse call
x[v,t] <- lsqnonlin(call_fit, x0, opts, y = rates[v,t], sigma = s[v])$x;
}
}
return(x)
}
################################################################################
# shadow rates, via Inverse Call Transformation
eta <- 0.005;
zeta <- 0;
icy <- InverseCallTransformation(y, tau, eta, zeta);
# log-rates
lny <- log(y)
# rate changes
dy <- diff(t(y));
dlny <- diff(t(log(y)));
dx <- diff(t(icy));
################################################################################
# plot rates/log-rates/shadow-rates #
################################################################################
fmt <- function(x) {
format(x, nsmall = 2, scientific = FALSE)
}
get_legend<-function(myggplot){
tmp <- ggplot_gtable(ggplot_build(myggplot))
leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
legend <- tmp$grobs[[leg]]
return(legend)
}
y_df <- melt(data.frame(date = as.Date(date), t(y)), id.vars = "date",
variable.name = "Maturity")
gp <- ggplot(y_df, aes(x = date, y = value)) +
geom_line(aes(colour = Maturity)) +
labs(list(title = "JGB Rates", x = "", y = "")) +
scale_x_date(expand = c(0, 0)) +
theme(legend.position = "bottom") +
theme(plot.margin = unit(c(0, .5, -0.5, .5), "lines"))
lny_df <- melt(data.frame(date = as.Date(date), t(lny)), id.vars = "date",
variable.name = "Maturity")
gp1 <- ggplot(lny_df, aes(x = date, y = value)) +
geom_line(aes(colour = Maturity)) +
labs(list(title = "JGB Log-Rates", x = "", y = "")) +
scale_x_date(expand = c(0, 0)) +
theme(legend.position="none") +
theme(plot.margin = unit(c(0, .5, -0.5, .5), "lines")) +
scale_y_continuous(labels = fmt)
icy_df <- melt(data.frame(date = as.Date(date), t(icy)), id.vars = "date",
variable.name = "Maturity")
gp2 <- ggplot(icy_df, aes(x = date, y = value)) +
geom_line(aes(colour = Maturity)) +
labs(list(title = "JGB Shadow Rates", x = "", y = "")) +
scale_x_date(expand = c(0, 0)) +
theme(legend.position = "none") +
theme(plot.margin = unit(c(0, .5, -0.5, .5), "lines"))
legend <- get_legend(gp)
gp <- gp + theme(legend.position="none")
grid.arrange(gp, gp1, gp2, legend, nrow=4, heights=c(3, 3, 3, 0.5)) 
The main point regarding the shadow rates is to…
Notice how the behavior of the unconstrained risk drivers attained by means of the inverse call transformation looks smoother than those of yields and log-yields, displaying a homogeneous evolution for both short and long-term series.
… which is further explained in Meucci’s “The Prayer”…
Homogeneity ensures that we can apply statistical techniques to the observed time series of the risk drivers {yt}t=1,…,T and project future distributions.
You can find out more about this quest for invariance in “The Prayer” at http://papers.ssrn.com/sol3/papers.cfm?abstract_id=1753788.
One question I’ll be sure to ask is how to model “real-life” negative interest rates that we observed recently in the Euro area – this Meucci model doesn’t appear to admit such…
If you are also taking the plunge and attending the course, do drop me an email and we’ll meet up!
Click here for the R code on GitHub.
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.