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.
Define number of contracts
curveLength = 12; % Add path to helper functions addpath Utilities
Load cleaned forward curve history data, futures returns statistics and implied volatilities.
load Data/FwdCurveHistory load Data/OptionVol
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));
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');
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
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 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