MATLAB to Mathematica CorrelationMatrix implementation

Question

This question focuses on translating/adapting custom defined MATLAB code for creating a CorrelationMatrix to Mathematica.

Software and hardware used:

• Mathematica Version 10.3.0.0
• Computer: Mac Pro (Late 2013)
• Processor: 3.5 GHz 6-Core Intel Xenon E5
• Memory: 32 GB 1866 MHzDDR3
• Graphics: AMD FirePro D500 3 GB

Background

I have begun the translation of a body of MATLAB code for Pricing Bermudan Swaptions in the LIBOR Market Model, from a dissertation by a Steffen Hippler.

From the paper's introduction:

The fundamental entity modelled in the LIBOR market model is the forward interest rate curve that specifies the simply compounded interest rate at any two given points in time.

...

Bermudan swaptions are interest rate derivatives with early exercise features that are among the most liquidly traded (exotic) interest rate derivative contracts. Consequently, their pricing and risk management is of high practical importance. The pricing of these instruments, however, poses significant conceptual and theoretical difficulties. This is due to the fact that, in general, the pricing in the LIBOR market model has to be carried out via Monte Carlo simulation techniques, which in turn do not lend themselves naturally to the pricing of options with early exercise features.

I hope to eventually post a demonstration of the implementation on the Wolfram Demonstrations Project site.

Note: If FinancialDerivate[] or CUDAFinancialDerivative[] enabled such pricing, it would greatly simply my objective.

It does not appear that either function supports Bermuda style options (options that one can only exercise at defined times) or swaptions, options on swaps of fixed for floating cashflows.

A possible solution to this larger objective (which would likely require specialized knowledge outside the general domain of this forum), could show a way to utilize existing FinancialDerivate[] or CUDAFinancialDerivative[] functionality to replicate pricing of Bermuda Swaptions.

As I translate and develop the larger solution, I want to take advantage of Mathematica's functional approach as well as its high-level functionality.

Background specific to the immediate question

Hippler includes the following two custom MatLab functions for building a correlation matrix:

function [corr] = CorrFunc(beta, horizon, i, j)
    M=horizon;
    corr=exp(-abs(i-j)/(M-1)*(-log(beta(3))+beta(1)*(i^2+j^2+i*j-3*M*i-3*M*j+3*i+3*j+2*M^2-M-4)*1/((M-2)*(M-3))-beta(2)*(i^2+j^2+i*j-M*i-M*j+3*i+3*j+3*M^2-2)*1/((M-2)*(M-3))));
end

function CORR = CorrelationMatrix(beta, maturities, horizon, startRate, endRate)
   CORR=zeros(endRate-startRate);
   for i = startRate:endRate
        for j = startRate:endRate
            CORR(i,j) = CorrFunc(beta, horizon, maturities(j), maturities(i));
        end 
   end
end

Some definitions/term-examples:

beta (β or beta coefficient) - of an investment provides a measure of an instrument's volatility relative to a specific market as a whole. An individual stock or bond may have more or less volatility than the the larger stock or bond markets.

horizon :

If you buy a zero-coupon bond and hold it to maturity, the ROR on your investment is the zero rate at which you bought the bond.

If you buy a t-year zero-coupon bond and sell it at time T

maturities - "maturity dates" on which the principal amount of the debt instruments become due and interest payments stop.

startRate, endRate - establishes a range of iteration.

The CorrFunc implements (in MATLAB) a parameterisation due to Schoenmakers and Coffey (J. Schoenmakers and C. Coffey, Systematic generation of parametric correlation structures for the libor market model, International Journal of Theoretical and Applied Finance (2002)) of the form:

$$\small \begin{align*}\rho_{ij}=\exp\left[-\frac{|i-j|}{M-1}\right.&\left(-\ln\beta_3+\beta_1\frac{i^2+j^2+ij-3Mi-3Mj+3i+3j+2M^2-M-4}{(M-2)(M-3)}\right.\\ &\left.\left.-\beta_2\frac{i^2+j^2+ij-Mi-Mj-3i-3j+3M^2-2}{(M-2)(M-3)}\right)\right]\end{align*}$$

with fitting parameter β = (β1, β2, β3) and M denoting the total number of rates under consideration.

I think, the MATLAB CorrFunc translates pretty directly to Mathematica:

corrFunc[beta_, horizon_, i_, j_] := Module[{M, corr},
  M = horizon;
  corr = Exp[-Abs[i - j]/(M - 1)*(-Log[beta[[3]]] + beta[[1]]*(i^2 + j^2 + i*j - 3*M*i - 3*M*j + 3*i + 3*j + 2*M^2 - M - 4)*1/((M - 2)*(M - 3)) - beta[[2]]*(i^2 + j^2 + i*j - M*i - M*j + 3*i + 3*j + 3*M^2 - 2)*1/((M - 2)*(M - 3)))]]

The MATLAB, CorrelationMatrix just looks like it simply applies the CorrFunc to a square matrix of the maturities.

So, it looks like I can simply do something like the following:

beta = {0.5, 0.4, 0.3};
horizon = 100;
Table[corrFunc[beta, horizon, i, j], {i, 10}, {j, 10}]

Note: I don't have MATLAB so I can't do a direct comparison of output.

Questions:

• Does it look like I'm on the right track?
• Whether, yes or no, can I alternatively do this by other/better means?
• LinearModelFit?
• GeneralizedLinearModelFit?
• Using Listable features?

As stated in the paper: The entire project of Pricing Bermudan Swaptions,

...poses significant conceptual and theoretical difficulties.

I want/hope to make use of all of Mathematica's advantages in simplifying code and writing efficient performing code.

It occurs to me that writing up the context of a question like this helps me, as an OP, better think about the problem. Not that I have a solution, but I hope it gets me into the area of the solution.

Answer 1

That was already a very good start. I take your function and modify it only slightly. I anticipate from the Matlab code that not the indices but the maturities should be arguments to corrFunc:

corrFunc[beta_, M_, mi_, mj_] := 
 Exp[-Abs[mi - mj]/(M - 1) (-Log[beta[[3]]] + 
     beta[[1]] (mi^2 + mj^2 + mi mj - 3 M mi - 3 M mj + 3 mi + 3 mj + 
        2 M^2 - M - 4) 1/((M - 2) (M - 3)) - 
     beta[[2]] (mi^2 + mj^2 + mi mj - M mi - M mj - 3 mi - 3 mj + 
        3 M^2 - 2) 1/((M - 2) (M - 3)))
  ]

Now we can loop directly over the list of maturities:

n = 100;
maturities = RandomReal[{-1, 1}, n];
beta = {0.5, 0.4, 0.3};
horizon = 100;
data0 = Table[
    corrFunc[beta, horizon, mi, mj], 
    {mi, maturities}, 
    {mj, maturities}
]; // AbsoluteTiming // First

0.181448

In order to speed that up, I created a CompiledFunction with similar argument pattern. For performance reasons, I substitute the list beta by three scalar arguments. Moreover, I apply N to the body of the function to convert all constants to reals such that the eventual C code is free of type casts.

ccorrFunc = Block[{mi, mj, M, β1, β2, β3},
 With[{code = N@corrFunc[{β1, β2, β3}, M, mi, mj]},
   Compile[
    {{β1, _Real}, {β2, _Real}, {β3, _Real}, {M, _Real}, {mi, _Real}, {mj, _Real}},
    code,
    CompilationTarget -> "C",
    RuntimeAttributes -> {Listable},
    Parallelization -> True,
    RuntimeOptions -> "Speed"
    ]
  ]
 ]

This is how we employ this new function to compute the matrix.

{mi, mj} = Transpose@ Flatten[Outer[List, maturities, maturities, 1], 1]; // AbsoluteTiming // First
data1 = Partition[
    ccorrFunc[Sequence @@ beta, horizon, mi, mj],
    n]; // AbsoluteTiming // First

0.00008

0.000518

This is about 300 times faster. And we obtain essentially the same result:

Max[Abs[data0 - data1]]

1.11022*10^-16

Since the input data is small compared to the actual computations, it might be also worthwhile to put that onto a graphics card. Unfortunately, my knowledge of OpenCL is a bit scarce... (CUDA would not work on the Fire Pro.)

Edit

Fixed signs in corrFun as suggested by J.M.

Answer 2

Starting from Henrik's implementation of corrFunc[], I have managed to shorten it slightly (after correcting a sign mistake):

corrFuncNu[beta_, M_, mi_, mj_] := 
Exp[-Abs[mi - mj] (Take[beta, 2].((mi^2 + mi mj + mj^2) {1, -1} +
                   {3 (1 - M), M + 3} (mi + mj) +
                   {2 M^2 - M - 4, 2 - 3 M^2})/((M - 2) (M - 3)) - Log[beta[[3]]])/(M - 1)]

We can then use DistanceMatrix[] with a custom distance function to generate the matrix:

n = 100;
BlockRandom[SeedRandom[42]; maturities = RandomReal[{-1, 1}, n]];
beta = {0.5, 0.4, 0.3}; horizon = 100;

DistanceMatrix[maturities, DistanceFunction -> Function[{mi, mj},
                           corrFuncNu[beta, horizon, mi, mj]]]

(in older versions, replace corrFuncNu[beta, horizon, mi, mj] with corrFuncNu[beta, horizon, First[mi], First[mj]])

On the computer I am using, this is twice as fast as generating the matrix with Table[], and just as accurate.