% Demand Forecast
% (Author Brian Kirtlan 26 July 2004)
% 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 excludes variation from forecast inputs

% SINGLE STAGE MODEL LOG VERSION
% DOMESTIC DEMAND IS TOTAL DEMAND BASED ON C + Residences + Price+ ShortageFlag (Price OK for log version as opposed to Total Energy cost)
% LIGHT IND AND COMM DEMAND IS TOTAL DEMAND BASED ON C + GDP + Shortage
% HEAVY DEMAND IS TOTAL DEMAND BASED ON TREND vs YEAR + Shortage (Non log)

% Model Standard Variables

MCRUNS = 100; %Number of Monte Carlo runs

CurrentCPIIndex = 1084; % 2003 figure (see input files for index)
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 2002) 
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

% 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)*CurrentCPIIndex;
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;

SingleRun = 0; % Set this to 1 to do a single model run using the raw input series (if 0 will run MonteCarlo)

% This section 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

if TrimTestOn == 1
    EndDom = EndDom - TrimPeriod;
    EndLightInd = EndLightInd - TrimPeriod;
    EndHeavyInd = EndHeavyInd - TrimPeriod; % (This test isn't very appropriate for Heavy Ind due to short data set)
    ForecastPeriods = TrimPeriod;
    TempFCGDP = FCGDP;
    TempFCHH = FCHH;
    TempFCLClosses = FCLClosses;
    TempFCNZPOP = FCNZPOP;
    TempFCShortage = FCShortage;
    TempFCYear = FCYear;
    TempFCPrice = FCPrice;
    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);
FittedDomDemand = zeros(EndDom-StartDom+1,1);
FittedDomDemandStore = zeros(EndDom-StartDom+1,MCRUNS);
FittedLightIndCommDemand = zeros(EndLightInd-StartLightInd+1,1);
FittedLightIndCommStore = zeros(EndLightInd-StartLightInd+1,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);

% temp storage for random inputs

GDPtemp = GDP;
TotalNZPoptemp = TotalNZPop;
DomResidencestemp = DomResidences;
DomRealPricetemp = DomRealPrice;

% Start of Monte Carlo

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

for MCCycle = 1:MCRUNS;

    % Create synthetic datasets (note - domestic residences assumed to vary in same
    % proportion as population)

    if SingleRun ~= 1;
        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));
            DomRealPrice(n)=DomRealPricetemp(n).*(1+normrnd(0,0.05)); % alternative of (1+PriceVariationProp(fix(rand(1)*(EndDom-PriceMovingAveragePeriod+1))+1)) too unstable;
        end;
    end;

    % Calculate Derived Series from input series

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


    % Estimate Domestic Model
    % THIS IS log(TOTAL DEMAND) BASED ON C + log(DomResidences) + log(Price) + Shortage

    DomesticModelData = [ones(size(DomDemand)) log(DomResidences) log(DomRealPrice) ShortageFlag];
    DomesticModelDataTrunc = DomesticModelData(StartDom:EndDom,:);
    DomDemTrunc = log(DomAdjDomDemand(StartDom:EndDom,:));
    DomDemandModel = DomesticModelDataTrunc\DomDemTrunc;

    % Domestic Fitted Results

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

    % Estimate Light Industrial and Commercial Model
    % THIS IS log(TOTAL DEMAND) BASED ON C + log(GDP) + Shortage

    LightIndCommModel = [ones(size(LightIndCommDemand)) log(GDP) ShortageFlag];
    LightIndCommDataTrunc = LightIndCommModel(StartLightInd:EndLightInd,:);
    LightIndCommDemandTrunc = log(LightIndCommDemand(StartLightInd:EndLightInd,:));
    LightIndCommModel = LightIndCommDataTrunc\LightIndCommDemandTrunc;

    % Light Industrial and Commercial Fitted Results

    FittedLightIndCommDemand = exp(LightIndCommDataTrunc*LightIndCommModel);
    FittedLightIndCommDemandComparison = [FittedLightIndCommDemand LightIndCommDemand(StartLightInd:EndLightInd,:)];

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

    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(FCHH(1:ForecastPeriods));
    DomDemandForecastInputs(:,3)=log(FCPrice(1:ForecastPeriods));
    DomDemandForecastInputs(:,4)=FCShortage(1:ForecastPeriods);
    DomTotalForecastDemand=exp(DomDemandForecastInputs*DomDemandModel);

    % Light Industrial and Commercial Forecast

    LightIndCommForecastInputs=zeros(ForecastPeriods,3);
    LightIndCommForecastInputs(:,1)=1;
    LightIndCommForecastInputs(:,2)=log(FCGDP(1:ForecastPeriods));
    LightIndCommForecastInputs(:,3)=FCShortage(1:ForecastPeriods);

    LightIndCommTotalForecastDemand = exp(LightIndCommForecastInputs*LightIndCommModel);

    % Heavy Industry Demand

    HeavyIndInputs = zeros(ForecastPeriods,3);
    HeavyIndInputs(:,1)=1;
    HeavyIndInputs(:,2)=FCYear;
    HeavyIndInputs(:,3)=FCShortage;
    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);
        [LightIndCommResultCheck, LightIndCommModelCheckCIs, LightIndCommErrorsCheck, LightIndCommErrorsCheckCIs, LightIndCommModelStats]=regress(LightIndCommDemandTrunc,LightIndCommDataTrunc);
        LightIndCommCoeffStats = tstats(LightIndCommDemandTrunc,LightIndCommDataTrunc);
        LightIndCommErrorsDW = DurbinWatson(LightIndCommErrorsCheck);

        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(LightIndCommModel)-1))./(EndLightInd-StartLightInd-(length(LightIndCommModel)-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;
    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;

    %Removed temporarily as getting errors due to the solutions that result in
    %rapid exponential growth
    %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);

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));
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:EndLightInd),FittedLightIndCommStore);
hold on;
LightIndCommDemandGraph=plot(Year(StartLightInd:EndLightInd),LightIndCommDemand(StartLightInd: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;
%hold off;

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

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;

% restore changed variables

GDP=GDPtemp;
TotalNZPop=TotalNZPoptemp;
DomResidences=DomResidencestemp;
DomRealPrice=DomRealPricetemp;

if TrimTestOn == 1
    FCGDP=TempFCGDP;
    FCHH = TempFCHH;
    FCLClosses = TempFCLClosses;
    FCNZPOP = TempFCNZPOP;
    FCShortage = TempFCShortage;
    FCYear = TempFCYear;
    FCPrice = TempFCPrice;
end;

%Warnings back on

warning on MATLAB:divideByZero
warning on MATLAB:nearlySingularMatrix

disp('Done')