spreadModel-dlm.R

#
# One-factor model of spread using DLM package
#
# Model in discrete form:
#
#   y[t] = a + b*s[t] + u[t]
#   s[t] = c*s[t-1] + v[t]
#
# Model in matrix form:
#
#                    _    _
#                   |      |
#   y[t] = [a, b] * |   1  | + u[t]
#                   |      |
#                   | s[t] |
#                   |_    _|
#
#    _    _       _    _     _      _    _     _
#   |      |     |      |   |        |   |      |
#   |   1  |     | 1  0 |   |   1    |   |   0  |
#   |      |  =  |      | * |        | + |      |
#   | s[t] |     | 0  c |   | s[t-1] |   | v[t] |
#   |_    _|     |_    _|   |_      _|   |_    _|
#
library(dlm)

#
# Given a set of DLM parameters, build a DLM model.
# The parameters are packed into a vector,
# so we must unpack all 7 element.
#
buildModelDLM = function(p) {
  a = p[[1]]
  b = p[[2]]
  c = p[[3]]
  dV = exp(p[[4]])
  dW = exp(p[[5]])
  s0 = p[[6]]
  sigma_s0 = exp(p[[7]])

  m0 = rbind(1, s0)
  C0 = matrix(c(0, 0,
                0, sigma_s0), 2, 2, byrow=T)

  FF = cbind(a, b)
  V = dV
  GG = matrix(c(1, 0,
                0, c), 2, 2, byrow=T)
  W = matrix(c(0, 0,
               0, dW), 2, 2, byrow=T)

  dlm(m0=m0, C0=C0, FF=FF, V=V, GG=GG, W=W)
}

#
# Given a time series, make a guess at the initial
# parameters for the DLM.
#
initialGuess = function(y) {
  a = 0; b = 1; c = 0.9
  dV = var(diff(y), na.rm=T)
  dW = dV / 2
  s0 = as.numeric(y[[1]])
  sigma_s0 = dV / 10

  c(a, b, c, log(dV), log(dW), s0, log(sigma_s0))
}

#
# Compute a typical spread:
# 10 yr Treasurys vs 5 yr Treasurys
#
library(xts)
library(quantmod)

yld10 = getSymbols("^TYX", auto.assign = F)
yld5 = getSymbols("^FVX", auto.assign = F)
y = Cl(yld10) - Cl(yld5)
colnames(y) = "Spread"

#
# Repeatedly call buildModelDLM() to find the MLE
#
mle = dlmMLE(y, parm=initialGuess(y), build=buildModelDLM)
if (mle$convergence != 0) stop(mle$message)

#
# From the MLE parameters, build the final model
#
model = buildModelDLM(mle$par)
cat("a =", model$FF[1,1], "\n")
cat("b =", model$FF[1,2], "\n")
cat("c =", model$GG[2,2], "\n")

#
# Plot the original data and the smoothed values
#
smooth.list = dlmSmooth(y, model)
smooth = xts::xts(smooth.list$s[-1,2], index(y))
plot.zoo(merge(y, smooth), screens=1,
         col=c("blue", "red"),
         main="Original and Smoothed Spread",
         xlab="Day", ylab="Spread (%)" )