# Natural Gas Storage Valuation: 3. Model Calibration and Validation

This is the third script of 4 in the Natural Gas Storage Valuation case study. In this script, we calibrate a forward curve model to historical and options-derived forward curve data. The forward curve model used is a multi-factor forward curve market model of the following form: dFn(t)/Fn(t) = mu(n) * dt + Sigma(n,:) * dZ where Fn is the n'th forward, mu is a vector of drifts, Sigma is a nForwards-by-nFactors matrix of volatility factors and dZ is a vector of nFactor uncorrelated brownian motions. Please see the Readme document for more information on the models that are used.

The model parameters include the drift *mu* and volatility factor matrix *Sigma*. These can be calibrated with both historical and forward-looking (option) information. The mu parameter is often set to 0 for risk-neutral valuation purposes.

In this script historical, as well as options-derived, means, covariances, volatilities and correlations are used to calibrate the model using principal components analysis to reduce the number of factors.

## Contents

## Model Structure

Define number of contracts

curveLength = 12; % Add path to helper functions addpath Utilities

## Import Forward Curve Data & Implied Volatilities

Load cleaned forward curve history data, futures returns statistics and implied volatilities.

load Data/FwdCurveHistory load Data/OptionVol

## Compute Curve Statistics

The forward curve model is a model of changes to the log of the forward curve. Compute these from the forward curve series, keeping only those maturities required for the model.

f = log(fwdPrice); df = diff(f, 1, 1); % First difference along first dimension df = collapseFwdMatrix(df); df = df(:, 1:curveLength); m = nanmean(df); sigma = nancov(df); vol = nanstd(df); cor = corr(df,'rows','pairwise'); visualizeVolatility(sigma, impVol(1:curveLength));

## Monthly Volatilities

Estimating a single set of volatilities & correlations for the aggregate set of futures returns averages away the seasonality in the data. We can instead compute monthly volatilities of returns grouped by the front month

[volm, corm, sigmam] = computeMonthlyVol(fwdPrice, expDate, 12); clf; bar3c(volm*sqrt(252)); ylabel('Front Month'); xlabel('Promptness'); zlabel('Volatility');

## Time Varying Correlations

Similar to the volatilities, the correlations between contracts (by promptness) change from month to month. This visual of monthly correlations demonstrates the need for a model that can account for this monthly variation

visualizeMonthlyCorr(corm, cor)

## Perform Principal Components Analysis

There is a lot of structure in the covariance matrix, suggesting that not all factors may be necessary to recreate the forward curve. Principal Components Analysis (PCA) can analyze this structure to produce a matrix of PCA coefficients and variances that specify the direction and magnitude of the most covariance in decending order. Usually only 2 or 3 principal components are necessary to explain most of the covariance.

nFactors = 3; [pcaCoeff, pcaVar, pctExpl] = pcacov(sigma); clf plot(cumsum(pctExpl), 'o-'); grid on; ylabel('Cumulative % Variance Explained'); xlabel('Number of Factors'); sigmaPCA = pcaCoeff(:,1:nFactors) * diag(sqrt(pcaVar(1:nFactors)));

Compare covariance matrix to the PCA-generated covariance matrix by selecting different number of factors

sigmaReco = sigmaPCA*sigmaPCA'; clf ax(1) = subplot(1,3,1, 'View', [103 28]); hold on; bar3c(sigma); axis tight; grid on; title('Original Covariance Matrix') ax(2) = subplot(1,3,2, 'View', [103 28]); hold on; bar3c(sigmaReco); axis tight; grid on title(sprintf('PCA Recreated Cov. Matrix with %d factors', nFactors)); ax(3) = subplot(1,3,3, 'View', [103 28]); hold on; err = sigma-sigmaReco; bar3c(err); axis tight; grid on; title('Error (Difference)'); linkprop(ax,'View'); fprintf('Maximum error between original and recovered sigma: %0.4f%%\n', max(max(abs(err)./sigma))*100);

Maximum error between original and recovered sigma: 3.7110%

## Save all calibrated models

Save calibrated parameters for use in subsequent analysis.

% Historical, Single sigma = calibratePCACov(sigma, 3); mu = m' + 1/2 * sum(sigma.^2, 2); save Models/FwdModel_HistPCA mu sigma % Forward Vols, Historical Corr, Single sigma = corr2cov(impVol(1:curveLength)/sqrt(252), cor); sigma = calibratePCACov(sigma, 3); mu = m' + 1/2 * sum(sigma.^2, 2); save Models/FwdModel_IVPCA mu sigma % Historical, Monthly sigma = calibratePCACov(sigmam, 3); mu = cellfun(@(sigma)m'+1/2*sum(sigma.^2, 2), sigma, 'UniformOutput', false); save Models/FwdModel_MonthlyPCA mu sigma