% New Zealand National Electricity Demand Forecast
% Author: Brian Kirtlan, New Zealand Electricity Commission
% Contact email: info@electricitycommission.govt.nz
% Phone : 0064 4 460 8860

% MonteCarlo Version - estimates model error by establishing a distribution around the long term trend for 
% the explanatory variables, and creating synthetic datasets which are used to re-estimate the model variables. 
% Forecasts are then calculated for each re-estimated model.

% This version includes variation from forecast input uncertainty 

% Model requires DemandForecastInputTablesDec2004.mat workspace

% Final Electricity Commission Proposed Model
% DOMESTIC DEMAND IS A SINGLE STAGE LOG MODEL: DEMAND/CAP BASED ON C + GDP/CAP + HH/CAP + Price
% LIGHT IND AND COMM DEMAND IS A TWO STAGE MODEL:  DEMAND BASED ON C + Year (1st stage only) + GDP + Shortage + lag
% HEAVY DEMAND IS A SINGLE STAGE MODEL:  TOTAL DEMAND BASED ON C + YEAR + Shortage

% Model Switches

SingleRun = 1; % Set this to 1 to do a single model run using the raw input series (if 0 will run MonteCarlo)
ForecastVariationOff = 0; %Set this to 1 to exclude forecast variation from Monte Carlo runs (if SingleRun is set to 1 then forecast variation is ignored anyway)
ModelVariationOff = 0; %Set this to 1 to exclude modelling uncertainty from the Monte Carlo runs (if SingleRun is set to 1 then this is ignored anway)

% These switches used for truncating the estimation period and forecasting across the
% remaining historical period - use the TrimTestOn var to switch the test on and off (it overwrites the normal forecasts) 

TrimTestOn = 0; %1 = run test, 0 = run normal forecast - i.e  this should normally be 0
TrimPeriod = 5; % Period to truncate models for when testing against historical data

% Model Standard Variables

MCRUNS = 1000; %Number of Monte Carlo runs

StartDom = 27; %Modelled Start Year of 1972 (less 1945 - the first year of data) 
EndDom = 58; %Historical data is current up to 2003 (less 1945)
StartLightInd = 26; %Modelled Start Year 1971 (less 1945)
EndLightInd = 58; %Data current to 2003 (less 1945)
StartHeavyInd = 44; %Modelled start year is 1989 (less 1945)
EndHeavyInd = 58; %Data current to 2003 (less 1945)
ForecastPeriods = 42; %Have forecast input data to 2045 (less 2003) 
ExistingEmbeddedGenerationGWh = 1337.1; % Gwh of currently available embedded gen - model assumes stays in proportion to total gen
StartHistDemand = 30; %First year that the combined totals of the historical demand data can be used for a total NZ historical timeseries is 1975 (less 1945)
EndHistDemand = 58; %Data current to 2003 (less 1945)
ForecastPercentile = 0.05; %Percentile at each tail used for forecast confidence limits (ie 0.05 = 90% confidence band)

%Spurious warnings off

warning off MATLAB:divideByZero
warning off MATLAB:nearlySingularMatrix

% temp storage for series used as random inputs or forecasts
% if the model falls over for some reason during a Monte Carlo run it means that the original data series in memory will need to be refreshed 

GDPtemp = GDP;
TotalNZPoptemp = TotalNZPop;
DomResidencestemp = DomResidences;
TempFCGDP = FCGDP;
TempFCHH = FCHH;
TempFCLClosses = FCLClosses;
TempFCNZPOP = FCNZPOP;
TempFCShortage = FCShortage;
TempFCYear = FCYear;
TempFCPrice = FCPrice;

% Establish input distributions

% GDP - calculate moving average and establish distribution table 
GDPMovingAveragePeriod = 5;
GDPMovingAverage=zeros((EndDom-GDPMovingAveragePeriod+1),1);
for i1 = 1:(EndDom-GDPMovingAveragePeriod+1);
    TempTotal = 0;
    for i2 = 0:(GDPMovingAveragePeriod-1);
        TempTotal = TempTotal+GDP(i1+i2);
    end;   
    GDPMovingAverage(i1)=TempTotal./GDPMovingAveragePeriod;
end;
GDPVariation=GDP((GDPMovingAveragePeriod-round(GDPMovingAveragePeriod/2)+1):(EndDom-round(GDPMovingAveragePeriod/2)+1))-GDPMovingAverage;
GDPVariationProp = GDPVariation./GDPMovingAverage;

% POP - calculate moving average and establish distribution table 
PopMovingAveragePeriod = 5;
PopMovingAverage=zeros((EndDom-PopMovingAveragePeriod+1),1);
for i1 = 1:(EndDom-PopMovingAveragePeriod+1);
    TempTotal = 0;
    for i2 = 0:(PopMovingAveragePeriod-1);
        TempTotal = TempTotal+TotalNZPop(i1+i2);
    end;   
    PopMovingAverage(i1)=TempTotal./PopMovingAveragePeriod;
end;
PopVariation=TotalNZPop((PopMovingAveragePeriod-round(PopMovingAveragePeriod/2)+1):(EndDom-round(PopMovingAveragePeriod/2)+1))-PopMovingAverage;
PopVariationProp = PopVariation./PopMovingAverage;

% Domestic Price - calculate moving average and establish distribution table
DomRealPrice=(DomElectricityPrice./CPI)*CPI(length(CPI));
DomRealPricetemp = DomRealPrice;
PriceMovingAveragePeriod=5;
PriceMovingAverage=zeros((EndDom-PriceMovingAveragePeriod+1),1);
for i1 = 1:(EndDom-PriceMovingAveragePeriod+1);
    TempTotal = 0;
    for i2 = 0:(PriceMovingAveragePeriod-1);
        TempTotal = TempTotal+DomRealPrice(i1+i2);
    end;
    PriceMovingAverage(i1)=TempTotal./PriceMovingAveragePeriod;
end;
PriceVariation=DomRealPrice((PriceMovingAveragePeriod-round(PriceMovingAveragePeriod/2)+1):(EndDom-round(PriceMovingAveragePeriod/2)+1))-PriceMovingAverage;
PriceVariationProp = PriceVariation./PriceMovingAverage;

if TrimTestOn == 1
    EndDom = EndDom - TrimPeriod;
    EndLightInd = EndLightInd - TrimPeriod;
    EndHeavyInd = EndHeavyInd - TrimPeriod; % (Note that this test isn't very appropriate for Heavy Ind due to short data set)
    ForecastPeriods = TrimPeriod;
    FCGDP = GDP((EndDom+1):(EndDom+TrimPeriod));
    FCHH = DomResidences((EndDom+1):(EndDom+TrimPeriod));
    FCLClosses = LineCoTransLosses((EndDom+1):(EndDom+TrimPeriod));
	FCNZPOP = TotalNZPop((EndDom+1):(EndDom+TrimPeriod));
	FCShortage = ShortageFlag((EndDom+1):(EndDom+TrimPeriod));
	FCYear = Year((EndDom+1):(EndDom+TrimPeriod));
	FCPrice = DomRealPrice((EndDom+1):(EndDom+TrimPeriod));
end;

if SingleRun == 1; MCRUNS = 1; end;

% Monte Carlo results stored in these

ForecastDemandStore = zeros(ForecastPeriods,MCRUNS); 
ForecastDomDemandStore = zeros(ForecastPeriods,MCRUNS);
ForecastLightIndCommDemandStore = zeros(ForecastPeriods,MCRUNS);
ForecastHeavyIndDemandStore = zeros(ForecastPeriods,MCRUNS);
ForecastGDPStore = zeros(ForecastPeriods,MCRUNS);
ForecastNZPOPStore = zeros(ForecastPeriods,MCRUNS);
ForecastHHStore = zeros(ForecastPeriods,MCRUNS);
ForecastPriceStore = zeros(ForecastPeriods,MCRUNS);
FittedDomDemandStore = zeros(EndDom-StartDom+1,MCRUNS);
FittedLightIndCommStore = zeros(EndLightInd-StartLightInd,MCRUNS); 
FittedHeavyIndStore = zeros(EndHeavyInd-StartHeavyInd+1,MCRUNS);
FittedTotalDemandStore = zeros(EndHeavyInd-StartHeavyInd+1,MCRUNS);
HistGDPStore = zeros(EndDom,MCRUNS);
HistPopStore = zeros(EndDom,MCRUNS);
HistPriceStore = zeros(EndDom,MCRUNS);
DomRsquaredStore = zeros(1,MCRUNS);
DomCorrectedRsquaredStore = zeros(1,MCRUNS);
LightIndCommRsquaredStore = zeros(1,MCRUNS);
LightIndCommCorrectedRsquaredStore = zeros(1,MCRUNS);
HeavyIndRsquaredStore = zeros(1,MCRUNS);
HeavyIndCorrectedRsquaredStore = zeros(1,MCRUNS);

% Plotting bin storage

NumberPlottingBins = 10;

TotalDemandBin = zeros(ForecastPeriods, NumberPlottingBins);
NormalisedBin = zeros(ForecastPeriods, NumberPlottingBins);

% Forecast input defined variation

GDPProductivityStdDev = 0.002; % Estimated based on NZIER assessment
FCNZPOPLastValueMean = 4920.09; % These are the mean and sigma values of a lognormal distribution fitted using the 2046 values of the
FCNZPOPLastValueLogSigma = 0.0785224; % various published Statistics New Zealand population scenarios (See PopForecasts.xls for data)

HouseSizeStdDev = 0.005; %Scenario based and phased in across total period over top of population driven changes
PriceStdDev = 0.1; % Estimated - ultimately this should be projected using the supply side modelling

% Calc forecast growth rate

FCGDPrate = FCGDP(2:ForecastPeriods)./FCGDP(1:ForecastPeriods-1);
FCNZPOPrate = FCNZPOP(2:ForecastPeriods)./FCNZPOP(1:ForecastPeriods-1);
FCHHrate = FCHH(2:ForecastPeriods)./FCHH(1:ForecastPeriods-1);

% Start of Monte Carlo

disp(['Starting Monte Carlo run with iterations set to ', int2str(MCRUNS)])

for MCCycle = 1:MCRUNS;

    % Create synthetic datasets

    if SingleRun ~= 1;
        % Select random adjustments for forecast variation that applies equally to all forecast years
        FCGDPProductivityVar = normrnd(0,GDPProductivityStdDev); %assumes normal distribution
        FCHouseSizeVar = normrnd(0,HouseSizeStdDev);
        FCPriceVar = normrnd(0,PriceStdDev);
        if (ModelVariationOff == 0);
            for n = 1:EndDom;
                GDP(n)=GDPtemp(n).*(1+GDPVariationProp(fix(rand(1)*(EndDom-GDPMovingAveragePeriod+1))+1));
                PopRand = rand(1);
                TotalNZPop(n)=TotalNZPoptemp(n).*(1+PopVariationProp(fix(PopRand*(EndDom-PopMovingAveragePeriod+1))+1));
                DomResidences(n)=DomResidencestemp(n).*(1+PopVariationProp(fix(PopRand*(EndDom-PopMovingAveragePeriod+1))+1)); %kept HH changes in proportion to pop changes 
                DomRealPrice(n)=DomRealPricetemp(n).*(1+normrnd(0,0.05)); % alternative of (1+PriceVariationProp(fix(rand(1)*(EndDom-PriceMovingAveragePeriod+1))+1)) too unstable;            
            end;
        end;
        FCEndPopNew = lognrnd(log(FCNZPOPLastValueMean), FCNZPOPLastValueLogSigma)*1000;
        FCNZPOPNew = TempFCNZPOP+(FCEndPopNew-TempFCNZPOP(ForecastPeriods))*(FCYear-FCYear(1)+1)./(FCYear(ForecastPeriods)-FCYear(1)+1);
        FCNZPOPnewrate = FCNZPOPNew(2:ForecastPeriods)./FCNZPOPNew(1:ForecastPeriods-1);
        FCNZPOPVar = FCNZPOPnewrate./FCNZPOPrate;
        if (ForecastVariationOff == 0) & (TrimTestOn~=1);
            % Calculate randomised forecast inputs
            for n = 2:ForecastPeriods;
                FCNZPOP(n) = FCNZPOP(n-1).*FCNZPOPnewrate(n-1);
                FCGDPnewrate = FCGDPrate(n-1).*(1+FCGDPProductivityVar).*(FCNZPOPVar(n-1)).*(1+GDPVariationProp(fix(rand(1)*(EndDom-GDPMovingAveragePeriod+1))+1));
                FCGDP(n)=FCGDP(n-1).*FCGDPnewrate;
                FCHHnewrate = FCHHrate(n-1).*(FCNZPOPVar(n-1));
                FCHH(n) = FCHH(n-1).*FCHHnewrate.*(1+(FCHouseSizeVar.*(n./ForecastPeriods))); %phases in the change in household size
                FCPrice(n) = TempFCPrice(n).*(1+FCPriceVar);
            end;        
        end;
    end;

    % Calculate Derived Series from input series

    DomGDPPerCap = GDP ./ TotalNZPop;
    DomAdjDomDemand = DomDemand + TempCorr;
    DomDemPerHH = DomAdjDomDemand ./ DomResidences .* 1000000;
    DomDemPerCap = DomAdjDomDemand ./ TotalNZPop;
    DomPeoplePerHH = TotalNZPop ./ DomResidences;
    DomHHPerCap = DomResidences ./ TotalNZPop;
    LightIndCommNomPrice = LightIndCommCost ./ LightIndCommDemand *100;
    LightIndCommRealPrice = LightIndCommNomPrice ./ CPI * CPI(length(CPI));
    LightIndCommDemandPerGDP = LightIndCommDemand ./ GDP;
    LightIndCommDemandPerCap = LightIndCommDemand ./ TotalNZPop;
    GDPPerCap = GDP ./ TotalNZPop;
    HeavyIndDemand = EDFIndCommDemand - LightIndCommDemand;

    % Estimate Domestic Model
    % THIS IS log(DEMAND PER PERSON) BASED ON C + log(GDP/CAP) and log(HH/CAP) + log(Price)

    DomesticModelData = [ones(size(DomDemand)) log(GDPPerCap) log(DomHHPerCap) log(DomRealPrice)];
    DomesticModelDataTrunc = DomesticModelData(StartDom:EndDom,:);
    DomDemTrunc = log(DomDemPerCap(StartDom:EndDom,:));
    DomDemandModel = DomesticModelDataTrunc\DomDemTrunc;

    % Domestic Fitted Results

    FittedDomDemand = exp((DomesticModelDataTrunc*DomDemandModel)).*TotalNZPop(StartDom:EndDom);
    FittedDomDemandComparison = [FittedDomDemand,DomAdjDomDemand(StartDom:EndDom)];

    % Estimate Light Industrial and Commercial Model
    % DEMAND  BASED ON C + Year (1st stage only) + GDP + Shortage + lag
    % Stage 1

    LightIndCommData = [ones(size(LightIndCommDemand)) Year GDP ShortageFlag];
    LightIndCommDataTrunc = LightIndCommData(StartLightInd:EndLightInd,:);
    Stage1LightIndCommDemand = LightIndCommDemand(StartLightInd:EndLightInd,:);
    Stage1LightIndCommModel = LightIndCommDataTrunc\Stage1LightIndCommDemand;
    Stage1LightIndCommResult = (LightIndCommDataTrunc*Stage1LightIndCommModel);

    % Stage 2

    LightIndCommStage2Trunc = [LightIndCommDataTrunc,Stage1LightIndCommResult];
    LightIndCommStage2Trunc(2:(EndLightInd-StartLightInd+1),5)=LightIndCommStage2Trunc(1:(EndLightInd-StartLightInd),5);
    LightIndCommStage2Trunc = LightIndCommStage2Trunc(2:(EndLightInd-StartLightInd+1),:);
    Stage2LightIndCommDemand = Stage1LightIndCommDemand(2:(EndLightInd-StartLightInd+1),:);
    LightIndCommStage2Trunc(:,2)=[];
    Stage2LightIndCommModel = LightIndCommStage2Trunc\Stage2LightIndCommDemand;

    % Light Industrial and Commercial Fitted Results

    FittedStage2LightIndCommResults = (LightIndCommStage2Trunc*Stage2LightIndCommModel);
    FittedLightIndCommDemand = FittedStage2LightIndCommResults;
    FittedLightIndCommDemandComparison = [FittedLightIndCommDemand LightIndCommDemand((StartLightInd+1):EndLightInd,:)];

    % Estimate Heavy Model
    % TOTAL DEMAND BASED ON TREND vs YEAR + ShortageFlag

    HeavyIndDataTrunc = [ones((EndHeavyInd-StartHeavyInd+1),1) Year(StartHeavyInd:EndHeavyInd) ShortageFlag(StartHeavyInd:EndHeavyInd)];
    HeavyIndDemandTrunc = HeavyIndDemand(StartHeavyInd:EndHeavyInd);
    HeavyIndModel = HeavyIndDataTrunc\HeavyIndDemandTrunc;
    HeavyIndResult = (HeavyIndDataTrunc*HeavyIndModel);
    FittedHeavyIndDemand = HeavyIndResult;
    FittedHeavyIndComparison = [FittedHeavyIndDemand HeavyIndDemandTrunc];

    % Forecast Demand Using Estimated Models

    % Domestic Forecast

    DomDemandForecastInputs=zeros(ForecastPeriods,4);
    DomDemandForecastInputs(:,1)=1;
    DomDemandForecastInputs(:,2)=log(FCGDP(1:ForecastPeriods)./FCNZPOP(1:ForecastPeriods));
    DomDemandForecastInputs(:,3)=log(FCHH(1:ForecastPeriods)./FCNZPOP(1:ForecastPeriods));
    DomDemandForecastInputs(:,4)=log(FCPrice(1:ForecastPeriods));
    DomTotalForecastDemand=exp((DomDemandForecastInputs*DomDemandModel)).*FCNZPOP(1:ForecastPeriods);

    % Light Industrial and Commercial Forecast

    LightIndCommForecastInputs = zeros(ForecastPeriods,4);
    LightIndCommForecastInputs(:,1) = 1;
    LightIndCommForecastInputs(:,2)=FCGDP(1:ForecastPeriods);
    LightIndCommForecastInputs(:,3)=FCShortage(1:ForecastPeriods);
    LightIndCommForecastInputs(1,4)=Stage1LightIndCommResult(EndLightInd-StartLightInd+1);
    LightIndCommForecast = zeros(ForecastPeriods,1);
    for i=1:ForecastPeriods;
        LightIndCommForecast(i)=LightIndCommForecastInputs(i,:)*Stage2LightIndCommModel;
        LightIndCommForecastInputs(i+1,4)=LightIndCommForecast(i);
    end;
    LightIndCommTotalForecastDemand = LightIndCommForecast;

    % Heavy Industry Demand

    HeavyIndInputs = zeros(ForecastPeriods,3);
    HeavyIndInputs(:,1)=1;
    HeavyIndInputs(:,2)=FCYear(1:ForecastPeriods);
    HeavyIndInputs(:,3)=FCShortage(1:ForecastPeriods);
    HeavyIndTotalForecastDemand = HeavyIndInputs*HeavyIndModel;

    % Forecast Lines Company Losses

    ForecastTotalLineLosses = FCLClosses.*(LightIndCommTotalForecastDemand + DomTotalForecastDemand);

    % Forecast Embedded Generation

    TotalCurrentDemand = (DomDemand(EndDom)+LightIndCommDemand(EndLightInd)+HeavyIndDemand(EndHeavyInd))+LineCoTransLosses(EndDom)*(DomDemand(EndDom)+LightIndCommDemand(EndLightInd));
    EmbeddedGenForecastDemand = ExistingEmbeddedGenerationGWh.*(HeavyIndTotalForecastDemand+LightIndCommTotalForecastDemand+DomTotalForecastDemand+ForecastTotalLineLosses)./TotalCurrentDemand;

    % Grand Total Forecast

    ForecastTotalNZDemand = DomTotalForecastDemand + LightIndCommTotalForecastDemand + HeavyIndTotalForecastDemand + ForecastTotalLineLosses - EmbeddedGenForecastDemand;

    % Estimated Historical GXP Demand

    HistTotalDemand = DomDemand(StartHistDemand:EndHistDemand) +LightIndCommDemand(StartHistDemand:EndHistDemand) + HeavyIndDemand(StartHistDemand:EndHistDemand);
    HistEstLocalLinesLosses = LineCoTransLosses(StartHistDemand:EndHistDemand) .* (DomDemand(StartHistDemand:EndHistDemand)+LightIndCommDemand(StartHistDemand:EndHistDemand));
    HistEstEmbedded = ExistingEmbeddedGenerationGWh*(HistTotalDemand+HistEstLocalLinesLosses)/(HistTotalDemand(EndHistDemand-StartHistDemand+1)+HistEstLocalLinesLosses(EndHistDemand-StartHistDemand+1));
    HistEstGXPDemand = HistTotalDemand + HistEstLocalLinesLosses - HistEstEmbedded;

    % Combine data for an actual and forecast graph

    CombinedDemand = [HistEstGXPDemand; ForecastTotalNZDemand];
    AllYears = [Year(StartHistDemand:EndHistDemand);FCYear];

    % Descriptive Stats (removing some during MonteCarlo as significantly slowing things down)

    if SingleRun == 1;
        DomDemandErrors = FittedDomDemandComparison(:,2)-FittedDomDemandComparison(:,1);
        [DomDemandResultsCheck, DomDemandResultsCheckCIs, DomDemandErrorsCheck, DomDemandErrorsCheckCIs, DomDemandModelStats] = regress(DomDemTrunc,DomesticModelDataTrunc);
        DomDemandCoeffStats = tstats((DomDemTrunc), DomesticModelDataTrunc);
        DomDemandErrorsDW = DurbinWatson(DomDemandErrorsCheck);

        LightIndCommErrors = FittedLightIndCommDemandComparison(:,2)-FittedLightIndCommDemandComparison(:,1);
        [Stage1LightIndCommResultCheck, Stage1LightIndCommModelCheckCIs, Stage1LightIndCommErrorsCheck, Stage1LightIndCommErrorsCheckCIs, Stage1LightIndCommModelStats]=regress(Stage1LightIndCommDemand,LightIndCommDataTrunc);
        Stage1LightIndCommCoeffStats = tstats(Stage1LightIndCommDemand,LightIndCommDataTrunc);
        Stage1LightIndCommErrorsDW = DurbinWatson(Stage1LightIndCommErrorsCheck);
        [Stage2LightIndCommResultCheck, Stage2LightIndCommModelCheckCIs, Stage2LightIndCommErrorsCheck, Stage2LightIndCommErrorsCheckCIs, Stage2LightIndCommModelStats]=regress(Stage2LightIndCommDemand,LightIndCommStage2Trunc);
        Stage2LightIndCommCoeffStats = tstats(Stage2LightIndCommDemand,LightIndCommStage2Trunc);
        Stage2LightIndCommErrorsDW = DurbinWatson(Stage2LightIndCommErrorsCheck);
        [Stage2LightIndCommErrorsDWh, Stage2LightIndCommErrorsDWhsig] = durbinh(Stage2LightIndCommErrorsDW,(Stage2LightIndCommCoeffStats(4,2))^2,length(Stage2LightIndCommErrorsCheck),1);

        HeavyIndErrors = FittedHeavyIndComparison(:,2)-FittedHeavyIndComparison(:,1);
        [HeavyIndResultsCheck,HeavyIndResultsCheckCIs,HeavyIndErrorsCheck,HeavyIndErrorsCheckCIs, HeavyIndModelStats] = regress(HeavyIndDemandTrunc,HeavyIndDataTrunc);
        HeavyIndCoeffStats = tstats(HeavyIndDemandTrunc,HeavyIndDataTrunc);
        HeavyIndErrorsDW = DurbinWatson(HeavyIndErrors);
    end;

    DomRSquaredFitted = FittedDomDemandComparison(:,1)-(sum(FittedDomDemandComparison(:,1)))/length(FittedDomDemandComparison(:,1));
    DomRSquaredActual = FittedDomDemandComparison(:,2)-(sum(FittedDomDemandComparison(:,2)))/length(FittedDomDemandComparison(:,2));
    DomDemandRsquared = (sum(DomRSquaredFitted.*DomRSquaredActual)./(((sum(DomRSquaredFitted.^2)).^0.5).*((sum(DomRSquaredActual.^2)).^0.5))).^2;
    DomDemandCorrectedRSquared = (((EndDom-StartDom)-1).*DomDemandRsquared-(length(DomDemandModel)-1))./(EndDom-StartDom-(length(DomDemandModel)-1)-1);

    LightIndCommRSquaredFitted = FittedLightIndCommDemandComparison(:,1)-(sum(FittedLightIndCommDemandComparison(:,1)))/length(FittedLightIndCommDemandComparison(:,1));
    LightIndCommRSquaredActual = FittedLightIndCommDemandComparison(:,2)-(sum(FittedLightIndCommDemandComparison(:,2)))/length(FittedLightIndCommDemandComparison(:,2));
    LightIndCommRsquared = (sum(LightIndCommRSquaredFitted.*LightIndCommRSquaredActual)./(((sum(LightIndCommRSquaredFitted.^2)).^0.5).*((sum(LightIndCommRSquaredActual.^2)).^0.5))).^2;
    LightIndCommCorrectedRSquared = (((EndLightInd-StartLightInd)-1).*LightIndCommRsquared-(length(Stage2LightIndCommModel)-1))./(EndLightInd-StartLightInd-(length(Stage2LightIndCommModel)-1)-1);

    HeavyIndRSquaredFitted = FittedHeavyIndComparison(:,1)-(sum(FittedHeavyIndComparison(:,1)))/length(FittedHeavyIndComparison(:,1));
    HeavyIndRSquaredActual = FittedHeavyIndComparison(:,2)-(sum(FittedHeavyIndComparison(:,2)))/length(FittedHeavyIndComparison(:,2));
    HeavyIndRsquared = (sum(HeavyIndRSquaredFitted.*HeavyIndRSquaredActual)./(((sum(HeavyIndRSquaredFitted.^2)).^0.5).*((sum(HeavyIndRSquaredActual.^2)).^0.5))).^2;
    HeavyIndCorrectedRSquared = (((EndHeavyInd-StartHeavyInd)-1).*HeavyIndRsquared-(length(HeavyIndModel)-1))./(EndHeavyInd-StartHeavyInd-(length(HeavyIndModel)-1)-1);

    % Store Monte Carlo data

    ForecastDemandStore(:,MCCycle)=ForecastTotalNZDemand;
    ForecastDomDemandStore(:,MCCycle)=DomTotalForecastDemand;
    ForecastLightIndCommDemandStore(:,MCCycle)=LightIndCommTotalForecastDemand;
    ForecastHeavyIndDemandStore(:,MCCycle)=HeavyIndTotalForecastDemand;
    ForecastGDPStore(:,MCCycle)=FCGDP;
    ForecastNZPOPStore(:,MCCycle)=FCNZPOP;
    ForecastHHStore(:,MCCycle)=FCHH;
    ForecastPriceStore(:,MCCycle)=FCPrice;
    FittedDomDemandStore(:,MCCycle)=FittedDomDemand;
    FittedLightIndCommStore(:,MCCycle)=FittedLightIndCommDemand;
    FittedHeavyIndStore(:,MCCycle)=FittedHeavyIndDemand;
    FittedTotalDemandStore(:,MCCycle)=FittedDomDemand((StartHeavyInd-StartDom):(EndDom-StartDom))+FittedLightIndCommDemand((StartHeavyInd-StartLightInd):(EndLightInd-StartLightInd))+FittedHeavyIndDemand;
    HistGDPStore(:,MCCycle)=GDP(1:EndDom);
    HistPopStore(:,MCCycle)=TotalNZPop(1:EndDom);
    HistPriceStore(:,MCCycle)=DomRealPrice(1:EndDom);
    DomRsquaredStore(1,MCCycle)=DomDemandRsquared;
    DomCorrectedRsquaredStore(1,MCCycle)=DomDemandCorrectedRSquared;
    LightIndCommRsquaredStore(1,MCCycle)=LightIndCommRsquared;
    LightIndCommCorrectedRsquaredStore(1,MCCycle)=LightIndCommCorrectedRSquared;
    HeavyIndRsquaredStore(1,MCCycle)=HeavyIndRsquared;
    HeavyIndCorrectedRsquaredStore(1,MCCycle)=HeavyIndCorrectedRSquared;

    for i=1:ForecastPeriods;
        TotalDemandBin(i,fix(ForecastTotalNZDemand(i)/(120000/NumberPlottingBins)))=TotalDemandBin(i,fix(ForecastTotalNZDemand(i)/(120000/NumberPlottingBins)))+1;
    end

end % of MonteCarlo

MaxTotalDemand = max(max(ForecastDemandStore)); 
maxcount = max(TotalDemandBin')';
for i=1:ForecastPeriods;
    NormalisedBin(i,:)=TotalDemandBin(i,:)./maxcount(i);
end

% Calculate forecast percentiles for each year

ForecastLowerPercentile = zeros(1,ForecastPeriods);
ForecastUpperPercentile = zeros(1,ForecastPeriods);
ForecastMedian = zeros(1,ForecastPeriods);
DomForecastLowerPercentile = zeros(1,ForecastPeriods);
DomForecastUpperPercentile = zeros(1,ForecastPeriods);
DomForecastMedian = zeros(1,ForecastPeriods);
LightIndCommForecastLowerPercentile = zeros(1,ForecastPeriods);
LightIndCommForecastUpperPercentile = zeros(1,ForecastPeriods);
LightIndCommForecastMedian = zeros(1,ForecastPeriods);
HeavyIndForecastLowerPercentile = zeros(1,ForecastPeriods);
HeavyIndForecastUpperPercentile = zeros(1,ForecastPeriods);
HeavyIndForecastMedian = zeros(1,ForecastPeriods);
NZPopForecastLowerPercentile = zeros(1,ForecastPeriods);
NZPopForecastUpperPercentile = zeros(1,ForecastPeriods);
NZPopForecastMedian = zeros(1,ForecastPeriods);
GDPForecastLowerPercentile = zeros(1,ForecastPeriods);
GDPForecastUpperPercentile = zeros(1,ForecastPeriods);
GDPForecastMedian = zeros(1,ForecastPeriods);
HHForecastLowerPercentile = zeros(1,ForecastPeriods);
HHForecastUpperPercentile = zeros(1,ForecastPeriods);
HHForecastMedian = zeros(1,ForecastPeriods);

if SingleRun~=1;
    SortedForecast = ForecastDemandStore;
    for i=1:ForecastPeriods;
        SortedForecast(i,:)=sort(SortedForecast(i,:));
    end;
    ForecastLowerPercentile = SortedForecast(:,fix(ForecastPercentile.*MCRUNS));
    ForecastUpperPercentile = SortedForecast(:,MCRUNS-fix(ForecastPercentile.*MCRUNS));
    ForecastMedian = SortedForecast(:,fix(MCRUNS./2));

    SortedForecast = ForecastDomDemandStore;
    for i=1:ForecastPeriods;
        SortedForecast(i,:)=sort(SortedForecast(i,:));
    end;
    DomForecastLowerPercentile = SortedForecast(:,fix(ForecastPercentile.*MCRUNS));
    DomForecastUpperPercentile = SortedForecast(:,MCRUNS-fix(ForecastPercentile.*MCRUNS));
    DomForecastMedian = SortedForecast(:,fix(MCRUNS./2));

    SortedForecast = ForecastLightIndCommDemandStore;
    for i=1:ForecastPeriods;
        SortedForecast(i,:)=sort(SortedForecast(i,:));
    end;
    LightIndCommForecastLowerPercentile = SortedForecast(:,fix(ForecastPercentile.*MCRUNS));
    LightIndCommForecastUpperPercentile = SortedForecast(:,MCRUNS-fix(ForecastPercentile.*MCRUNS));
    LightIndCommForecastMedian = SortedForecast(:,fix(MCRUNS./2));

    SortedForecast = ForecastHeavyIndDemandStore;
    for i=1:ForecastPeriods;
        SortedForecast(i,:)=sort(SortedForecast(i,:));
    end;
    HeavyIndForecastLowerPercentile = SortedForecast(:,fix(ForecastPercentile.*MCRUNS));
    HeavyIndForecastUpperPercentile = SortedForecast(:,MCRUNS-fix(ForecastPercentile.*MCRUNS));
    HeavyIndForecastMedian = SortedForecast(:,fix(MCRUNS./2));

    SortedForecast = ForecastNZPOPStore;
    for i=1:ForecastPeriods;
        SortedForecast(i,:)=sort(SortedForecast(i,:));
    end;
    NZPopForecastLowerPercentile = SortedForecast(:,fix(ForecastPercentile.*MCRUNS));
    NZPopForecastUpperPercentile = SortedForecast(:,MCRUNS-fix(ForecastPercentile.*MCRUNS));
    NZPopForecastMedian = SortedForecast(:,fix(MCRUNS./2));

    SortedForecast = ForecastGDPStore;
    for i=1:ForecastPeriods;
        SortedForecast(i,:)=sort(SortedForecast(i,:));
    end;
    GDPForecastLowerPercentile = SortedForecast(:,fix(ForecastPercentile.*MCRUNS));
    GDPForecastUpperPercentile = SortedForecast(:,MCRUNS-fix(ForecastPercentile.*MCRUNS));
    GDPForecastMedian = SortedForecast(:,fix(MCRUNS./2));

    SortedForecast = ForecastHHStore;
    for i=1:ForecastPeriods;
        SortedForecast(i,:)=sort(SortedForecast(i,:));
    end;
    HHForecastLowerPercentile = SortedForecast(:,fix(ForecastPercentile.*MCRUNS));
    HHForecastUpperPercentile = SortedForecast(:,MCRUNS-fix(ForecastPercentile.*MCRUNS));
    HHForecastMedian = SortedForecast(:,fix(MCRUNS./2));
end;

if SingleRun == 1;
    ColumnNames = char('Var      ','       Coeff','    SD','    T Value','  P');    
    DomRowNames = char('Constant           ', 'log(GDPPerCap)', 'log(DomHHPerCap)', 'log(DomRealPrice)');
    LightIndCommRowNames = char('Constant           ','GDP','Shortage','Lag');
    HeavyIndRowNames = char('Constant           ','Year','Shortage');
    disp([' ']);
    disp(['Domestic Model Results']);
    disp([' ']);
    disp(cellstr(ColumnNames)');
    disp([DomRowNames num2str(DomDemandCoeffStats,5)]);
    disp([' ']);
    disp(['Light Industrial and Commercial Model Results']);
    disp([' ']);
    disp(cellstr(ColumnNames)');
    disp([LightIndCommRowNames num2str(Stage2LightIndCommCoeffStats,5)]);
    disp([' ']);
    disp(['Heavy Industrial Model Results']);
    disp([' ']);
    disp(cellstr(ColumnNames)');
    disp([HeavyIndRowNames num2str(HeavyIndCoeffStats,5)]);
end

% Graph results 

figure(1);
FCgraph=plot(FCYear(1:ForecastPeriods),ForecastDemandStore);
axis([FCYear(1) FCYear(ForecastPeriods) 0 200000]);
hold on;
title('Forecast Total Demand');
UpperPCGraph = plot(FCYear(1:ForecastPeriods),ForecastUpperPercentile,'k-');
LowerPCGraph = plot(FCYear(1:ForecastPeriods),ForecastLowerPercentile,'k-');
MedianForecastGraph = plot(FCYear(1:ForecastPeriods),ForecastMedian,'k-');
set(UpperPCGraph,'LineWidth',3);
set(LowerPCGraph,'LineWidth',3);
set(MedianForecastGraph,'LineWidth',3);
drawnow;
hold off;

%figure(2);
%GDPgraph=plot(Year(1:EndDom),HistGDPStore);
%hold on;
%GDPActualGraph=plot(Year(1:EndDom),GDPtemp,'k-');
%set(GDPActualGraph,'LineWidth',2);
%title('GDP synthetic input series');
%drawnow;
%hold off;

%figure(3);
%Popgraph=plot(Year(1:EndDom),HistPopStore);
%hold on;
%PopActualGraph=plot(Year(1:EndDom),TotalNZPoptemp,'k-');
%set(PopActualGraph,'LineWidth',2);
%title('Population synthetic input series');
%drawnow;
%hold off;

figure(4);
DomFitgraph=plot(Year(StartDom:EndDom),FittedDomDemandStore);
hold on;
DomDemandGraph=plot(Year(StartDom:EndDom),DomDemand(StartDom:EndDom),'k-');
set(DomDemandGraph,'LineWidth',2);
title('Fitted Domestic Demand');
drawnow;
hold off;

figure(5);
LightIndFitgraph=plot(Year(StartLightInd+1:EndLightInd),FittedLightIndCommStore);
hold on;
LightIndCommDemandGraph=plot(Year(StartLightInd+1:EndLightInd),LightIndCommDemand(StartLightInd+1:EndLightInd),'k-');
set(LightIndCommDemandGraph,'LineWidth',2);
title('Fitted Light Industrial and Commercial Demand');
drawnow;
hold off;

figure(6);
HeavyIndgraph=plot(Year(StartHeavyInd:EndHeavyInd),FittedHeavyIndStore);
hold on;
HeavyIndDemandGraph=plot(Year(StartHeavyInd:EndHeavyInd),HeavyIndDemand(StartHeavyInd:EndHeavyInd),'k-');
set(HeavyIndDemandGraph,'LineWidth',2);
title('Fitted Heavy Industrial Demand');
drawnow;

%figure(7);
%TotalFittedgraph=plot(FittedTotalDemandStore);
%hold on;
%TotalDemandGraph=plot(HistTotalDemand(15:(EndHistDemand-StartHistDemand+1)),'k-');
%set(TotalDemandGraph,'LineWidth',2);
%title('Fitted Total Demand');
%drawnow;

figure(8);
hist(DomRsquaredStore);
title('Domestic Model Rsquared');
drawnow;

figure(9);
hist(LightIndCommRsquaredStore);
title('Light Industrial and Commercial Rsquared');
drawnow;

figure(10);
hist(HeavyIndRsquaredStore);
title('Heavy Industrial Rsquared');
drawnow;

figure(11);
DomFCgraph=plot(FCYear(1:ForecastPeriods),ForecastDomDemandStore);
if TrimTestOn == 1; MaxY = 15000;  else MaxY = 30000; end; 
axis([FCYear(1) FCYear(ForecastPeriods) 0 MaxY]);
hold on;
title('Forecast Domestic Demand');
UpperDomPCGraph = plot(FCYear(1:ForecastPeriods),DomForecastUpperPercentile,'k-');
LowerDomPCGraph = plot(FCYear(1:ForecastPeriods),DomForecastLowerPercentile,'k-');
MedianDomForecastGraph = plot(FCYear(1:ForecastPeriods),DomForecastMedian,'k-');
set(UpperDomPCGraph,'LineWidth',3);
set(LowerDomPCGraph,'LineWidth',3);
set(MedianDomForecastGraph,'LineWidth',3);
if TrimTestOn ==1
    ActualDomTrim = plot(FCYear(1:ForecastPeriods),DomDemand((EndDom+1):(EndDom+TrimPeriod)),'r-');
    set(ActualDomTrim,'LineWidth',3);
end;
drawnow;
hold off;

figure(12);
LightIndCommFCgraph=plot(FCYear(1:ForecastPeriods),ForecastLightIndCommDemandStore);
if TrimTestOn == 1; MaxY = 20000;  else MaxY = 50000; end; 
axis([FCYear(1) FCYear(ForecastPeriods) 0 MaxY]);
hold on;
title('Forecast Light Industrial and Commercial Demand');
UpperLightIndCommPCGraph = plot(FCYear(1:ForecastPeriods),LightIndCommForecastUpperPercentile,'k-');
LowerLightIndCommPCGraph = plot(FCYear(1:ForecastPeriods),LightIndCommForecastLowerPercentile,'k-');
MedianLightIndCommForecastGraph = plot(FCYear(1:ForecastPeriods),LightIndCommForecastMedian,'k-');
set(UpperLightIndCommPCGraph,'LineWidth',3);
set(LowerLightIndCommPCGraph,'LineWidth',3);
set(MedianLightIndCommForecastGraph,'LineWidth',3);
if TrimTestOn ==1
    ActualLightIndTrim = plot(FCYear(1:ForecastPeriods),LightIndCommDemand((EndLightInd+1):(EndLightInd+TrimPeriod)),'r-');
    set(ActualLightIndTrim,'LineWidth',3);
end;
drawnow;
hold off;

figure(13);
HeavyIndFCgraph=plot(FCYear(1:ForecastPeriods),ForecastHeavyIndDemandStore);
if TrimTestOn == 1; MaxY = 10000;  else MaxY = 20000; end;
axis([FCYear(1) FCYear(ForecastPeriods) 0 MaxY]);
hold on;
title('Forecast Heavy Industrial Demand');
UpperHeavyIndPCGraph = plot(FCYear(1:ForecastPeriods),HeavyIndForecastUpperPercentile,'k-');
LowerHeavyIndPCGraph = plot(FCYear(1:ForecastPeriods),HeavyIndForecastLowerPercentile,'k-');
MedianHeavyIndForecastGraph = plot(FCYear(1:ForecastPeriods),HeavyIndForecastMedian,'k-');
set(UpperHeavyIndPCGraph,'LineWidth',3);
set(LowerHeavyIndPCGraph,'LineWidth',3);
set(MedianHeavyIndForecastGraph,'LineWidth',3);
if TrimTestOn ==1
    ActualHeavyIndTrim = plot(FCYear(1:ForecastPeriods),HeavyIndDemand((EndHeavyInd+1):(EndHeavyInd+TrimPeriod)),'r-');
    set(ActualHeavyIndTrim,'LineWidth',3);
end;
drawnow;
hold off;

% The graph produced in figure 14 has embedded gen as a negative so that the total demand is correct. 
% The embedded gen is then repasted over the graph a second time otherwise it is overwritten by residential
% The down side of this is that the residential looks smaller in the graph than it actually is.

figure(14);
ForecastTotalNZDemandComponents = [(EmbeddedGenForecastDemand.*-1) DomTotalForecastDemand LightIndCommTotalForecastDemand HeavyIndTotalForecastDemand ForecastTotalLineLosses];
TotalDemandGraph = area(FCYear, ForecastTotalNZDemandComponents);
axis([FCYear(1) FCYear(ForecastPeriods) -10000 100000]);
legend('Embedded Gen','Residential','Comm & Light Ind','Heavy Ind','Line Losses');
hold on;
title('Forecast Total New Zealand Demand');
TotalDemandGraph = area(FCYear,(EmbeddedGenForecastDemand.*-1));
hold off;
drawnow;

figure(15);
ForecastGDPGraph = plot(FCYear,ForecastGDPStore);
axis([FCYear(1) FCYear(ForecastPeriods) 0 300000]);
hold on;
title('Forecast GDP');
UpperGDPPCGraph = plot(FCYear(1:ForecastPeriods),GDPForecastUpperPercentile,'k-');
LowerGDPPCGraph = plot(FCYear(1:ForecastPeriods),GDPForecastLowerPercentile,'k-');
MedianGDPForecastGraph = plot(FCYear(1:ForecastPeriods),GDPForecastMedian,'k-');
set(UpperGDPPCGraph,'LineWidth',3);
set(LowerGDPPCGraph,'LineWidth',3);
set(MedianGDPForecastGraph,'LineWidth',3);
if TrimTestOn ==1
    ActualGDPTrim = plot(FCYear(1:ForecastPeriods),GDP((EndDom+1):(EndDom+TrimPeriod)),'r-');
    set(ActualGDPTrim,'LineWidth',3);
end;
drawnow;
hold off;

figure(16);
ForecastPOPGraph = plot(FCYear,ForecastNZPOPStore);
axis([FCYear(1) FCYear(ForecastPeriods) 0 7000000]);
title('Forecast Population');
hold on;
UpperPopPCGraph = plot(FCYear(1:ForecastPeriods),NZPopForecastUpperPercentile,'k-');
LowerPopPCGraph = plot(FCYear(1:ForecastPeriods),NZPopForecastLowerPercentile,'k-');
MedianPopForecastGraph = plot(FCYear(1:ForecastPeriods),NZPopForecastMedian,'k-');
set(UpperPopPCGraph,'LineWidth',3);
set(LowerPopPCGraph,'LineWidth',3);
set(MedianPopForecastGraph,'LineWidth',3);
if TrimTestOn ==1
    ActualPopTrim = plot(FCYear(1:ForecastPeriods),TotalNZPop((EndDom+1):(EndDom+TrimPeriod)),'r-');
    set(ActualPopTrim,'LineWidth',3);
end;
drawnow;
hold off;

figure(17);
ForecastHHGraph = plot(FCYear,ForecastHHStore);
axis([FCYear(1) FCYear(ForecastPeriods) 0 3000000]);
title('Forecast Households');
drawnow;

figure(18);
ForecastPriceGraph = plot(FCYear,ForecastPriceStore);
axis([FCYear(1) FCYear(ForecastPeriods) 0 25]);
title('Forecast Residential Price');
drawnow;

%figure(19);
%ForecastGDPGrowthRate = ForecastGDPStore(2:ForecastPeriods,:)./ForecastGDPStore(1:ForecastPeriods-1,:)-1;
%ForecastGDPGrowthRateGraph = plot(FCYear(2:ForecastPeriods), ForecastGDPGrowthRate);
%axis([FCYear(1) FCYear(ForecastPeriods) 0 0.05]);
%drawnow;

figure(20);
FCgraph2=plot(FCYear(1:ForecastPeriods),ForecastDemandStore);
axis([Year(StartHistDemand) FCYear(ForecastPeriods) 0 120000]);
hold on;
HistDemandgraph = plot(Year(StartHistDemand:EndHistDemand), HistEstGXPDemand,'k-');
title('Forecast Total Demand');
if SingleRun ~= 1
    HistToForecastGap = plot([Year(EndHistDemand); FCYear(1)], [HistEstGXPDemand(EndHistDemand-StartHistDemand+1); ForecastMedian(1)], 'k-');
    UpperPCGraph = plot(FCYear(1:ForecastPeriods),ForecastUpperPercentile,'k-');
    LowerPCGraph = plot(FCYear(1:ForecastPeriods),ForecastLowerPercentile,'k-');
    MedianForecastGraph = plot(FCYear(1:ForecastPeriods),ForecastMedian,'k-');
    set(UpperPCGraph,'LineWidth',3);
    set(LowerPCGraph,'LineWidth',3);
    set(MedianForecastGraph,'LineWidth',3);
    set(HistToForecastGap,'LineWidth',3);
end;
set(HistDemandgraph,'LineWidth',3);
set(gca, 'Ygrid','on');
drawnow;
hold off;

% restore changed variables

GDP = GDPtemp;
TotalNZPop = TotalNZPoptemp;
DomResidences = DomResidencestemp;
FCGDP = TempFCGDP;
FCHH = TempFCHH;
FCLClosses = TempFCLClosses;
FCNZPOP = TempFCNZPOP;
FCShortage = TempFCShortage;
FCYear = TempFCYear;
FCPrice = TempFCPrice;

%Warnings back on

warning on MATLAB:divideByZero
warning on MATLAB:nearlySingularMatrix

disp('Done')