MultiFactorForwardCurveSimulator(nSims, Seed, historicalForwardCurves, nDays, nFactors)

function [FCSims, AllVolFNs, volFNsCorrs, initFCs]  = MultiFactorForwardCurveSimulator(nSims, Seed, historicalForwardCurves, nDays, nFactors)

%**************************************************************************
%
% This module simulates the Clewlow and Strickland multi factor 
% forward curves for multi commodities.
%
% The calibration is automated by using the full historical dataset.
%
% Inputs:
%        nSims - number of simulations
%        Seed - seed for random number generator
%        FC - Contains historical forward data
%        nDays - number of days to simulate
%        nFactors - number of factors to use
%
% Outputs:
%        FCSims - daily FC sims
%        AllVolFNs - all the volatility functions
%        volFNsCorrs - volatility functions
%        initFCs - last FCs to start the simulations
%
% Ahmos Sansom - March 2013
%**************************************************************************

% Determine number of markets, historical prices and maturities
NMarkets = size(historicalForwardCurves,3);
NMaturities =  size(historicalForwardCurves,2);
NPrices =  size(historicalForwardCurves,1);


%% Step 1 - Calibration

% Allocate
returns = zeros(NPrices -1, NMaturities, NMarkets);
CoVEachMarket = zeros(NMaturities, NMaturities, NMarkets);
AllVolFNs = zeros(NMaturities, NMaturities, NMarkets);

for i = 1:NMarkets
    % 1.1 Calculate relative log returns
    returns(:,:, i) = price2ret(historicalForwardCurves(:,:,i),1:NPrices,'Continuous');
    
    % sometimes returns are normalised to zero mean - see discussion in
    % Carol Alexander and http://web.mit.edu/ceepr/www/publications/workingpapers/2009-002.pdf
    % Ignored here as well are seasonally affects and returns are relative,
    % i.e. aligned to front month - improvemnts can be to ensure that when
    % the contract changes, the returns are aligned accordingly, could add
    % exponential weighting...etc

    % 1.2 - PCA of each commodity
    % get covariance matrix of returns
    CoVEachMarket(:,:,i) = cov(returns(:,:,i));

    % PCA of covariance matrix of each commodity
    [coeff, latent] = pcacov(CoVEachMarket(:,:,i));
    
    %Determine volatility functions
    volFn = zeros(NMaturities, NMaturities);
    for j=1:NMaturities
        for k=1:NMaturities
            volFn(k,j) = coeff(k,j) * sqrt(latent(j,1));
        end
    end
    AllVolFNs(:,:,i) = volFn;
end

% 1.3 Determine correlation matrix across volatility functions

% The idea is to mininmise [V]_{A}[\rho][V]_{B}^{T}-[\Sigma]_{AB} = 0 to solve
% for the correlation matrix [\rho].
%
% where: [V]_{A} is the set of volatility function for market A and
% [V]_{B} is the set of volatility function for market B.
%
% The initial correlation matrix is based on the initial correlation across the factors.
% 
% See Multi-Factor Mulit-Commodity Models & Parameter Estimation Process by Lacima March 2008 that can improve the initial starting point and 
% http://web.mit.edu/ceepr/www/publications/workingpapers/2009-002.pdf for some technical explanation.

%Allocate
volFNsCorrs = zeros(NMarkets*nFactors);
tmpCorrelationMatrix = zeros(nFactors);

%Optimisation parameters - upper and lower bounds
options = optimset('TolFun',1e-10,'TolX',1e-10);
lb = -ones(nFactors);
ub = ones(nFactors);

for i=1:NMarkets
    for j=i:NMarkets
        for k=1:nFactors
            for l=1:nFactors
                if(i==j)
                    % volatility functions in the same market, so 0 or 1 (diagonal). 
                    if (k==l)
                       volFNsCorrs(nFactors*(i-1) + k, nFactors*(j-1) + l) = 1;
                    else
                       volFNsCorrs(nFactors*(i-1) + k, nFactors*(j-1) + l) = 0;
                    end
                else
                    tmpCorrelationMatrix(k, l) = corr(AllVolFNs(:,k,i),AllVolFNs(:,l,j));
                    % Determine correlation and take account of symmetry
                    volFNsCorrs(nFactors * (i-1) + k, nFactors * (j-1) + l) = tmpCorrelationMatrix(k, l);
                    volFNsCorrs(nFactors * (j-1) + l, nFactors * (i-1) + k) = volFNsCorrs(nFactors * (i-1) + k, nFactors * (j-1) + l);
                end
            end
        end
        
        % Optimise correlation matrix                
        if (i~=j )
            covAB = zeros(NMaturities);
            for k=1:NMaturities
                for l=1:NMaturities
                    covAB(k,l) = GetCov(returns(:,k, i),returns(:,l, j));
                    
                end
            end
            %covAB = returns(:,:, i)' * returns(:,:, j); % if mean is zero
            optimisedCorrs = lsqnonlin(@(x)f(x,AllVolFNs(:,1:nFactors,i),AllVolFNs(:,1:nFactors,j)',covAB),tmpCorrelationMatrix,lb,ub,options);
            
            % add the optimised correlation matrix
            for k=1:nFactors
                for l=1:nFactors
                    volFNsCorrs(nFactors * (i-1) + k, nFactors * (j-1) + l) = optimisedCorrs(k, l);
                    volFNsCorrs(nFactors * (j-1) + l, nFactors * (i-1) + k) = volFNsCorrs(nFactors * (i-1) + k, nFactors * (j-1) + l);
                end
            end
        end
                        
        
    end
end

% Probably should add some code to ensure correlation matrix is positive
% semi definite.

%% Step 2 - Simulate

% Check number of days to simulate is valid - assume 21 working days per
% matuirity
if(nDays > NMaturities*21)
  error(message('Too may days for forward curve'));
end

% initialise random number generator
s = RandStream('twister','Seed',Seed);
RandStream.setGlobalStream(s)

% Get last forward curves - starting point for simulations
initFCs = zeros(NMarkets,NMaturities);
for i=1:NMarkets
    initFCs(i,:) = historicalForwardCurves(end,:,i);
end

%Allocate simulated forward sims variable
FCSims = zeros(nSims,NMaturities,nDays,NMarkets);

% Set up matrix with number of working days per maturity - assume that each
% month has 21 working days per maturity
workingDaysperMaturity = zeros(NMaturities,1);
count = 0;
for i=1:NMaturities
    workingDaysperMaturity(i) = 21;
    count = count + workingDaysperMaturity(i);
    if(count > nDays)
        workingDaysperMaturity(i) = nDays -(i-1)*21;
        break;
    end
end


% Loop plan
        % 1. First two loops control the number of days to time march the FC; NMaturities is the number of maturitis in the FC
        %    2. Each maturity has n days to maturity which is also a loop
        %       3. Loop over the number of simulations (nSims)
        %          4. Loop over every maturity (iMaturity always starting from first month to ensure that each day moving forward, the volatility structured is rolled forward.
        %             5.  Loop over every forward curve market(NMarkets)
        %                6. Loop over evey factor (nFactors) to determine volatility and drfit
        %                

% Loop round the maturity steps, i.e how many days to run for each maturity
monthStartDay = 1;
for iMonth = 1:NMaturities
            
    %Number of days to simulate, i.e. number of days until contract matures for each month                
    ndaysPerMat = workingDaysperMaturity(iMonth);
               
    % j is the number days to simulate over, starts from 1
    for j = monthStartDay:(ndaysPerMat + monthStartDay - 1)
        j % output day
        
        %Get correlated random numbers using cholesky decomposition
        mydeltaz = randn(nSims,nFactors*NMarkets)*chol(volFNsCorrs);
        
        for i = 1:nSims % simulate curves
                                        
            % Loop round maturities - reduces as you move into the next
            % maturity bucket
            for iMaturity = 1:NMaturities - (iMonth - 1)
                                                    
                 % Set maturity step to store results in correct location -
                 % assumes move into next maturity using iMonth, could
                 % assume iMatStep = iMaturity, but forward curve would
                 % always start from same point as moving forward.
                 iMatStep = iMaturity + (iMonth - 1);

                 for iCurve = 1:NMarkets % Loop around markets
                    
                     %Initialise volatility 
                     sigma1 = 0.0;
                     sigma2 = 0.0;

                     for k = 1:nFactors
                         % Assumes volatility function always start from
                         % first vol (iMaturity) when moving forward in
                         % time to next maturity bucket and time step is
                         % set to unity
                         sigma1 = sigma1 + AllVolFNs(iMaturity, k, iCurve) * AllVolFNs(iMaturity, k, iCurve);
                         sigma2 = sigma2 + AllVolFNs(iMaturity, k, iCurve) * mydeltaz(i, nFactors * (iCurve - 1) + k);
                     end                                                        % Pick up respective correlated random number
                                
                     % update forward curve
                     if (j == 1)
                         FCSims(i, iMatStep, j, iCurve) = initFCs(iCurve, iMatStep) * exp(-0.5 * sigma1 + sigma2);
                     else
                         FCSims(i, iMatStep, j, iCurve) = FCSims(i, iMatStep, j - 1, iCurve) * exp(-0.5 * sigma1 + sigma2);
                     end
                     
                 end
            end
        end % end of sims
    end
    
    %updated counter
    monthStartDay = monthStartDay + ndaysPerMat;
    
end