# 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.

## Model Structure

Define number of contracts

```curveLength = 12;
% Add path to helper functions
```

## Import Forward Curve Data & Implied Volatilities

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

```load Data/FwdCurveHistory
```

## 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)');

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
```