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

% Last amended December 2009

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

% Model requires DemandForecastInputTablesSept2007.mat workspace to be in same directory as the script

% DOMESTIC DEMAND IS A LOG MODEL: DEMAND/CAP BASED ON C + GDP/CAP + HH/CAP + Price
% INDUSTRIAL AND COMMERCIAL DEMAND IS LINEAR MODEL:  DEMAND BASED ON C + GDP + Shortage
% HEAVY INDUSTRIAL (TIWAI)  : TIWAI DEMAND FORECAST AS CONSTANT (MAXIMUM OF HISTORICAL ANNUAL DEMAND)

clear all

load('2009EnergyInputsHistoric.mat');
load('2009EnergyInputsForecast.mat');

OutputDir = 'C:\DataFiles\MATLABWORK\NatDemand2009\FinalFiles\NatDemand2009\EnergyModel\';
OutputFN = 'DFResults.xls';


% Model Switches

SingleRun = 1; % Set this to 1 to do a single model run using the raw input series (if 0 will run MonteCarlo)
MCRUNS = 100; % Number of Monte Carlo runs
ForecastVariationOff = 0; % Set this to 1 to exclude forecast variation from the Monte Carlo runs (if SingleRun is set to 1 then forecast variation is ignored)
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)
ShowGraphs = 0; % Set to 1 to show graphs

% 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

StartYrResidentialMod = 1974; % First year to use historical data from for estimating the residential model 
StartYrIndCommMod = 1986; % First year to use historical data from for estimating the industrial/commercial model 
StartYrHeavyIndMod = 1974; % First year to use historical data from for estimating the heavy industrial model 
DataYrEnd = Year(length(Year)); % Last year that historical data covers unless overwritten
MaxForecastYear = FCYear(length(FCYear)); % Last year that forecast data covers unless overwritten

ExistingEmbeddedGenerationGWh = constants(6); % Gwh of currently available embedded gen (excludes cogen so that consistent with MED consumption figures) - model assumes forecast total embedded stays in proportion to total demand
ForecastPercentile = constants(7); % Percentile at each tail used for forecast confidence limits (ie 0.05 = 90% confidence band, 0.1 = 80% confidence band)
PropDirectConnect = constants(8); % This is the current proportion of non-Tiwai load that is direct connect (used to estimate future local lines losses) - See regional model for source of figure - all the precision is to try to keep the regional and national model national results as close as possible

% Find first input data and forecast data year

FirstHistYear = Year(1);
FirstFCYear = FCYear(1);

StartDom = StartYrResidentialMod-FirstHistYear+1; % Array position of starting year based on selected start year 
EndDom = DataYrEnd-FirstHistYear+1; % Array position of last year of historical data to use in model
StartIndComm = StartYrIndCommMod-FirstHistYear+1; % Array position of starting year based on selected start year 
EndIndComm = DataYrEnd-FirstHistYear+1; % Array position of last year of historical data to use in model
StartHeavyInd = StartYrHeavyIndMod-FirstHistYear+1; % array position of starting year based on selected start year 
EndHeavyInd = DataYrEnd-FirstHistYear+1; % Array position of last year of historical data to use in model
ForecastPeriods = MaxForecastYear-FirstFCYear+1; % Array position of last year to forecast to
StartHistDemand = 1974-FirstHistYear+1; % First year to use in total NZ historical timeseries (data currently runs across all categoires from 1974)
EndHistDemand = DataYrEnd-FirstHistYear+1; % Array position of last year of historical data used in model

HistFittedStart = max([StartDom, StartIndComm, StartHeavyInd]);% sets position for the first year that total fitted historical demand can be calculated for

% Model start date error catching

if StartDom < 1; error('Residential model start date earlier than historical input data'); end;
if StartIndComm < 1; error('Industrial/commerical model start date earlier than historical input data'); end;
if StartHeavyInd < 1; error('Heavy industrial model start date earlier than historical input data'); end;
if StartHistDemand < 1; error('Date specified for combined historical series earlier than historical input data'); end;

% Forecast input defined variation

GDPProductivityStdDev = constants(1); % Estimated based on NZIER assessment
FCNZPOPLastValueLogMean = constants(2); % These are the mean and sigma values of a lognormal distribution fitted using the 2050 values of the
FCNZPOPLastValueLogSigma = constants(3);  % various published Statistics New Zealand population scenarios (See Population projections summary.xls for data)
HouseSizeStdDev = constants(4); %Scenario based and phased in across total period over top of population driven changes
PriceStdDev = constants(5); % Estimated - ultimately this should be projected using the supply side modelling

tic

% Spurious warnings off - causes spam during some Monte Carlo or trucation runs depending on period selection 

warning off MATLAB:divideByZero
warning off MATLAB:nearlySingularMatrix

% 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 

% temp storage for series used as random inputs or forecasts
GDPtemp = GDP;
TotalNZPoptemp = TotalNZPop;
DomResidencestemp = DomResidences;
TempFCGDP = FCGDP;
TempFCHH = FCHH;
TempFCLClosses = FCLClosses;
TempFCNZPOP = FCNZPOP;
TempFCShortage = FCShortage;
TempFCYear = FCYear;
TempFCPrice = FCPrice;

% Seed the random number generator
state = 123;
rand('state', state);
randn('state', state);

% 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;
GDPVariationPropInit = GDPVariation./GDPMovingAverage;
GDPVariationProp = GDPVariationPropInit-(sum(GDPVariationPropInit)/length(GDPVariationPropInit)); % the two tails of the moving average affect the mean of this distribution so it needs to be adjusted so that total variation had mean of zero 

% 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;
PopVariationPropInit = PopVariation./PopMovingAverage;
PopVariationProp = PopVariationPropInit-(sum(PopVariationPropInit)/length(PopVariationPropInit));

% 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;
PriceVariationPropInit = PriceVariation./PriceMovingAverage;
PriceVariationProp = PriceVariationPropInit-(sum(PriceVariationPropInit)/length(PriceVariationPropInit));

if TrimTestOn == 1
    EndDom = EndDom - TrimPeriod;
    EndIndComm = EndIndComm - TrimPeriod;   
    EndHeavyInd = EndHeavyInd - TrimPeriod; % (In this case this just sets the final year for the range that the max load is calculated off for forecasting forward)
    EndHistDemand = EndHistDemand - TrimPeriod;
    DataYrEnd = DataYrEnd - TrimPeriod;
    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);
ForecastIndCommDemandStore = 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);
FittedIndCommStore = zeros(EndIndComm-StartIndComm+1,MCRUNS);
FittedHeavyIndStore = zeros(EndHeavyInd-StartHeavyInd+1,MCRUNS);
FittedTotalDemandStore = zeros(DataYrEnd-FirstHistYear-HistFittedStart+2,MCRUNS);
HistGDPStore = zeros(EndDom,MCRUNS);
HistPopStore = zeros(EndDom,MCRUNS);
HistPriceStore = zeros(EndDom,MCRUNS);
DomRsquaredStore = zeros(1,MCRUNS);
DomCorrectedRsquaredStore = zeros(1,MCRUNS);
IndCommRsquaredStore = zeros(1,MCRUNS);
IndCommCorrectedRsquaredStore = zeros(1,MCRUNS);
HeavyIndRsquaredStore = zeros(1,MCRUNS);
HeavyIndCorrectedRsquaredStore = zeros(1,MCRUNS);

% Plotting bin storage

NumberPlottingBins = 200;

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

% 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(' ');
disp('National demand forecast model');
disp(' ');
disp('*******************');
disp(['Starting Monte Carlo run with iterations set to ', int2str(MCRUNS)]);

for MCCycle = 1:MCRUNS;

    if mod(MCCycle,1000)==0; disp(['Run   ' num2str(MCCycle)]); toc; end; % display when hit every 1000 runs
    
    % 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(FCNZPOPLastValueLogMean, FCNZPOPLastValueLogSigma);
        FCNZPOPNew = TempFCNZPOP(1:ForecastPeriods)+(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*(n./ForecastPeriods)));
            end;        
        end;
    end;

    % Calculate Derived Series from input series

    DomGDPPerCap = GDP ./ TotalNZPop;
    DomDemPerHH = DomDemand ./ DomResidences .* 1000000;
    DomDemPerCap = DomDemand ./ TotalNZPop;
    DomPeoplePerHH = TotalNZPop ./ DomResidences;
    DomHHPerCap = DomResidences ./ TotalNZPop;
       
    IndCommDemand = EDFIndCommDemand - TiwaiDemand;
    HeavyIndDemand = TiwaiDemand;
    
    GDPPerCap = GDP ./ TotalNZPop;
    
    % 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,DomDemand(StartDom:EndDom)];

    % Estimate Industrial and Commercial Model
    % DEMAND  BASED ON C + GDP + Shortage

    IndCommData = [ones(size(IndCommDemand)) GDP ShortageFlag];
    IndCommDataTrunc = IndCommData(StartIndComm:EndIndComm,:);
    IndCommDemandTrunc = IndCommDemand(StartIndComm:EndIndComm,:);
    IndCommModel = IndCommDataTrunc\IndCommDemandTrunc;
    
    % Industrial and Commercial Fitted Results

    FittedIndCommDemand = IndCommDataTrunc*IndCommModel;
    FittedIndCommDemandComparison = [FittedIndCommDemand IndCommDemandTrunc];
   
    % Heavy Industrial Model
    % TIWAI TOTAL DEMAND CONSTANT BASED ON MAXIMUM DEMAND 
    
    HeavyIndDemandTrunc = HeavyIndDemand(StartHeavyInd:EndHeavyInd);
    HeavyIndConstant = max(HeavyIndDemandTrunc);
    FittedHeavyIndDemand = HeavyIndDemandTrunc; % Set fit to actual
    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);

    % Industrial and Commercial Forecast

    IndCommForecastInputs = zeros(ForecastPeriods,3);
    IndCommForecastInputs(:,1) = 1;
    IndCommForecastInputs(:,2) = FCGDP(1:ForecastPeriods);
    IndCommForecastInputs(:,3) = FCShortage(1:ForecastPeriods);
    IndCommTotalForecastDemand = IndCommForecastInputs*IndCommModel;
    
    % Heavy Industrial Demand
         
    HeavyIndTotalForecastDemand (1:ForecastPeriods,1) = HeavyIndConstant;
    
    % Forecast Lines Company Losses, Embedded Gen and Total Demand
  
    ForecastTotalLineLosses = FCLClosses.*(IndCommTotalForecastDemand + DomTotalForecastDemand).*(1-PropDirectConnect);
    TotalCurrentDemand = (DomDemand(EndDom)+IndCommDemand(EndIndComm)+HeavyIndDemand(EndHeavyInd))+LineCoTransLosses(EndDom)*(DomDemand(EndDom)+IndCommDemand(EndIndComm));    
    EmbeddedGenForecastDemand = ExistingEmbeddedGenerationGWh.*(HeavyIndTotalForecastDemand+IndCommTotalForecastDemand+DomTotalForecastDemand+ForecastTotalLineLosses)./TotalCurrentDemand;
    ForecastTotalNZDemand = DomTotalForecastDemand + IndCommTotalForecastDemand + HeavyIndTotalForecastDemand + ForecastTotalLineLosses - EmbeddedGenForecastDemand;
        
    % Estimated Historical GXP Demand

    HistTotalDemand = DomDemand(StartHistDemand:EndHistDemand) + IndCommDemand(StartHistDemand:EndHistDemand) + HeavyIndDemand(StartHistDemand:EndHistDemand);
    HistEstLocalLinesLosses = LineCoTransLosses(StartHistDemand:EndHistDemand) .* (DomDemand(StartHistDemand:EndHistDemand)+IndCommDemand(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);

        IndCommErrors = FittedIndCommDemandComparison(:,2)-FittedIndCommDemandComparison(:,1);
        [IndCommResultCheck, IndCommModelCheckCIs, IndCommErrorsCheck, IndCommErrorsCheckCIs, IndCommModelStats]=regress(IndCommDemandTrunc,IndCommDataTrunc);
        IndCommCoeffStats = tstats(IndCommDemandTrunc,IndCommDataTrunc);
        IndCommErrorsDW = DurbinWatson(IndCommErrorsCheck);
    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);

    IndCommRSquaredFitted = FittedIndCommDemandComparison(:,1)-(sum(FittedIndCommDemandComparison(:,1)))/length(FittedIndCommDemandComparison(:,1));
    IndCommRSquaredActual = FittedIndCommDemandComparison(:,2)-(sum(FittedIndCommDemandComparison(:,2)))/length(FittedIndCommDemandComparison(:,2));
    IndCommRSquared = (sum(IndCommRSquaredFitted.*IndCommRSquaredActual)./(((sum(IndCommRSquaredFitted.^2)).^0.5).*((sum(IndCommRSquaredActual.^2)).^0.5))).^2;
    IndCommCorrectedRSquared = (((EndIndComm-StartIndComm)-1).*IndCommRSquared-(length(IndCommModel)-1))./(EndIndComm-StartIndComm-(length(IndCommModel)-1)-1);
    
    % Store Monte Carlo data

    ForecastDemandStore(:,MCCycle)=ForecastTotalNZDemand;
    ForecastDomDemandStore(:,MCCycle)=DomTotalForecastDemand;
    ForecastIndCommDemandStore(:,MCCycle)=IndCommTotalForecastDemand;
    ForecastHeavyIndDemandStore(:,MCCycle)=HeavyIndTotalForecastDemand;
    ForecastGDPStore(:,MCCycle)=FCGDP;
    ForecastNZPOPStore(:,MCCycle)=FCNZPOP;
    ForecastHHStore(:,MCCycle)=FCHH;
    ForecastPriceStore(:,MCCycle)=FCPrice;
    FittedDomDemandStore(:,MCCycle)=FittedDomDemand;
    FittedIndCommStore(:,MCCycle)=FittedIndCommDemand;
    FittedHeavyIndStore(:,MCCycle)=FittedHeavyIndDemand;
    FittedTotalDemandStore(:,MCCycle)=FittedDomDemand((HistFittedStart-StartDom+1):(EndDom-StartDom+1))+FittedIndCommDemand((HistFittedStart-StartIndComm+1):(EndIndComm-StartIndComm+1))+HistEstLocalLinesLosses(HistFittedStart-StartHistDemand+1:EndHistDemand-StartHistDemand+1)-HistEstEmbedded(HistFittedStart-StartHistDemand+1:EndHistDemand-StartHistDemand+1)+FittedHeavyIndDemand((HistFittedStart-StartHeavyInd+1):(EndHeavyInd-StartHeavyInd+1));
    HistGDPStore(:,MCCycle)=GDP(1:EndDom);
    HistPopStore(:,MCCycle)=TotalNZPop(1:EndDom);
    HistPriceStore(:,MCCycle)=DomRealPrice(1:EndDom);
    DomRsquaredStore(1,MCCycle)=DomDemandRSquared;
    DomCorrectedRsquaredStore(1,MCCycle)=DomDemandCorrectedRSquared;
    IndCommRsquaredStore(1,MCCycle)=IndCommRSquared;
    IndCommCorrectedRsquaredStore(1,MCCycle)=IndCommCorrectedRSquared;
    
    for i=1:ForecastPeriods;
        TotalDemandBin(i,fix(ForecastTotalNZDemand(i)/(150000/NumberPlottingBins)))=TotalDemandBin(i,fix(ForecastTotalNZDemand(i)/(150000/NumberPlottingBins)))+1;
    end

end % of MonteCarlo

% xlswrite code
ExcelOutput1 = [FCYear,DomTotalForecastDemand, IndCommTotalForecastDemand,HeavyIndTotalForecastDemand,ForecastTotalLineLosses,EmbeddedGenForecastDemand,ForecastTotalNZDemand];
xlswrite(strcat(OutputDir,OutputFN),ExcelOutput1,'NatForecast', 'a10');


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(ForecastPeriods,1);
ForecastUpperPercentile = zeros(ForecastPeriods,1);
ForecastMedian = zeros(ForecastPeriods,1);
DomForecastLowerPercentile = zeros(ForecastPeriods,1);
DomForecastUpperPercentile = zeros(ForecastPeriods,1);
DomForecastMedian = zeros(ForecastPeriods,1);
IndCommForecastLowerPercentile = zeros(ForecastPeriods,1);
IndCommForecastUpperPercentile = zeros(ForecastPeriods,1);
IndCommForecastMedian = zeros(ForecastPeriods,1);
HeavyIndForecastLowerPercentile = zeros(ForecastPeriods,1);
HeavyIndForecastUpperPercentile = zeros(ForecastPeriods,1);
HeavyIndForecastMedian = zeros(ForecastPeriods,1);
NZPopForecastLowerPercentile = zeros(ForecastPeriods,1);
NZPopForecastUpperPercentile = zeros(ForecastPeriods,1);
NZPopForecastMedian = zeros(ForecastPeriods,1);
GDPForecastLowerPercentile = zeros(ForecastPeriods,1);
GDPForecastUpperPercentile = zeros(ForecastPeriods,1);
GDPForecastMedian = zeros(ForecastPeriods,1);
HHForecastLowerPercentile = zeros(ForecastPeriods,1);
HHForecastUpperPercentile = zeros(ForecastPeriods,1);
HHForecastMedian = zeros(ForecastPeriods,1);

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 = ForecastIndCommDemandStore;
    for i=1:ForecastPeriods;
        SortedForecast(i,:)=sort(SortedForecast(i,:));
    end;
    IndCommForecastLowerPercentile = SortedForecast(:,fix(ForecastPercentile.*MCRUNS));
    IndCommForecastUpperPercentile = SortedForecast(:,MCRUNS-fix(ForecastPercentile.*MCRUNS));
    IndCommForecastMedian = 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)');
    IndCommRowNames = char('Constant           ','GDP','Shortage');
    CommRowNames = char('Constant           ','GDP','Shortage');
    IndRowNames = char('Constant           ','GDP','Shortage');
    HeavyIndRowNames = char('Constant           ','GDP','Shortage');
    disp(' ');
    disp('Domestic Model Results');
    disp(' ');
    disp(cellstr(ColumnNames)');
    disp([DomRowNames num2str(DomDemandCoeffStats,5)]);
    disp(['R-squared   ' num2str(DomDemandRSquared)]);
    disp(['R-squared Adj  ' num2str(DomDemandCorrectedRSquared)]);
    disp(['DurbinWatson  ' num2str(DomDemandErrorsDW)]);
    disp(' ');   
    disp('Industrial and Commercial Model Results');
    disp(' ');
    disp(cellstr(ColumnNames)');
    disp([IndCommRowNames num2str(IndCommCoeffStats,5)]);
    disp(['R-squared   ' num2str(IndCommRSquared)]);
    disp(['R-squared Adj  ' num2str(IndCommCorrectedRSquared)]);
    disp(['DurbinWatson  ' num2str(IndCommErrorsDW)]);
end

% Graph results 

if ShowGraphs == 1
    
    % Index of figures
    disp(' ');
    disp ('Index of displayed graphs');
    disp(' ');
    disp(' Figure 1 : Forecast total demand');
    disp(' Figure 2 : GDP input series');
    disp(' Figure 3 : Population input series');
    disp(' Figure 4 : Fitted residential demand');
    disp(' Figure 5 : Fitted industrial and commercial demand');
    disp(' Figure 6 : Fitted heavy industrial demand');
    disp(' Figure 7 : Fitted total demand');
    disp(' Figure 8 : Domestic model r-squared');
    disp(' Figure 9 : Industrial and commercial r-squared');
    disp(' Figure 10 : Heavy industrial r-squared');
    disp(' Figure 11 : Forecast domestic demand');
    disp(' Figure 12 : Forecast industrial and commercial demand');
    disp(' Figure 13 : Forecast heavy industrial demand');
    disp(' Figure 14 : Forecast total demand by component');
    disp(' Figure 15 : Forecast GDP');
    disp(' Figure 16 : Forecast population');
    disp(' Figure 17 : Forecast households');
    disp(' Figure 18 : Forecast residential price');
    disp(' Figure 19 : Forecast GDP % annual growth');
    disp(' Figure 20 : Forecast total demand with historical demand');
    disp(' Figure 21 : Forecast total demand with historical demand by component');
    disp(' Figure 22 : Forecast total demand by component with embedded and line losses combined');
    
    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(1:EndDom),'k-');
    set(GDPActualGraph,'LineWidth',2);
    title('GDP input series');
    drawnow;
    hold off;

    figure(3);
    Popgraph=plot(Year(1:EndDom),HistPopStore);
    hold on;
    PopActualGraph=plot(Year(1:EndDom),TotalNZPoptemp(1:EndDom),'k-');
    set(PopActualGraph,'LineWidth',2);
    title('Population 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);
    IndCommFitgraph=plot(Year(StartIndComm:EndIndComm),FittedIndCommStore);
    hold on;
    IndCommDemandGraph=plot(Year(StartIndComm:EndIndComm),IndCommDemand(StartIndComm:EndIndComm),'k-');
    set(IndCommDemandGraph,'LineWidth',2);
    title('Fitted 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(Year(HistFittedStart:EndHistDemand),FittedTotalDemandStore);
    hold on;
    TotalDemandGraph=plot(Year(HistFittedStart:EndHistDemand),HistTotalDemand(HistFittedStart-StartHistDemand+1:(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(IndCommRsquaredStore);
    title('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);
    IndCommFCgraph=plot(FCYear(1:ForecastPeriods),ForecastIndCommDemandStore);
    if TrimTestOn == 1; MaxY = 22000;  else MaxY = 60000; end; 
    axis([FCYear(1) FCYear(ForecastPeriods) 0 MaxY]);
    hold on;
    title('Forecast Industrial and Commercial Demand');
    UpperIndCommPCGraph = plot(FCYear(1:ForecastPeriods),IndCommForecastUpperPercentile,'k-');
    LowerIndCommPCGraph = plot(FCYear(1:ForecastPeriods),IndCommForecastLowerPercentile,'k-');
    MedianIndCommForecastGraph = plot(FCYear(1:ForecastPeriods),IndCommForecastMedian,'k-');
    set(UpperIndCommPCGraph,'LineWidth',3);
    set(LowerIndCommPCGraph,'LineWidth',3);
    set(MedianIndCommForecastGraph,'LineWidth',3);
    if TrimTestOn ==1
        ActualIndCommTrim = plot(FCYear(1:ForecastPeriods),IndCommDemand((EndIndComm+1):(EndIndComm+TrimPeriod)),'r-');
        set(ActualIndCommTrim,'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);
    ForecastTotalNZDemandComponents1 = [(EmbeddedGenForecastDemand.*-1) DomTotalForecastDemand IndCommTotalForecastDemand HeavyIndTotalForecastDemand ForecastTotalLineLosses];
    TotalDemandGraph = area(FCYear, ForecastTotalNZDemandComponents1);
    axis([FCYear(1) FCYear(ForecastPeriods) -10000 100000]);
    legend('Embedded Gen','Residential','Ind & Comm','Heavy Ind','Distribution Line Losses');
    hold on;
    title('Forecast Total New Zealand Demand');
    ylabel('GWh');
    TotalDemandGraph = area(FCYear,(EmbeddedGenForecastDemand.*-1));
    hold off;
    drawnow;

    figure(15);
    ForecastGDPGraph = plot(FCYear,ForecastGDPStore);
    axis([FCYear(1) FCYear(ForecastPeriods) 0 400000]);
    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);
    title('Forecast GDP annual % growth');
    axis([FCYear(1) FCYear(ForecastPeriods) -0.05 0.1]);
    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;

    figure(21);
    if TrimTestOn == 0;    
        StartCombinedGraphYear = 1993;
        StartCombOffset = StartCombinedGraphYear-Year(StartHistDemand)+1;
        EndCombHist = EndHistDemand-StartHistDemand+1;
        CombinedTotalNZDemandComponents = [HistEstEmbedded(StartCombOffset:EndCombHist).*-1 DomDemand(StartHistDemand+StartCombOffset-1:EndHistDemand) IndCommDemand(StartHistDemand+StartCombOffset-1:EndHistDemand) HeavyIndDemand(StartHistDemand+StartCombOffset-1:EndHistDemand) HistEstLocalLinesLosses(StartCombOffset:EndCombHist);...
            (EmbeddedGenForecastDemand.*-1) DomTotalForecastDemand IndCommTotalForecastDemand HeavyIndTotalForecastDemand ForecastTotalLineLosses];
        CombinedYears = [(Year(StartHistDemand+StartCombOffset-1):Year(EndHistDemand))';FCYear];
        TotalDemandGraph = area(CombinedYears, CombinedTotalNZDemandComponents);
        axis([CombinedYears(1) CombinedYears(length(CombinedYears)) -10000 100000]);
        legend('Embedded Gen','Residential','Comm & Light Ind','Heavy Ind','Distribution Line Losses');
        hold on;
        title('Forecast Total New Zealand Demand');
        ylabel('GWh');
        TotalDemandGraph = area(CombinedYears,[HistEstEmbedded(StartCombOffset:EndCombHist).*-1; EmbeddedGenForecastDemand.*-1]);
        drawnow;
        hold off;
    end;

    figure(22);
    ForecastOther =   ForecastTotalLineLosses - EmbeddedGenForecastDemand;
    ForecastTotalNZDemandComponents = [DomTotalForecastDemand IndCommTotalForecastDemand HeavyIndTotalForecastDemand ForecastOther];
    TotalDemandGraph = area(FCYear, ForecastTotalNZDemandComponents);
    axis([FCYear(1) FCYear(ForecastPeriods) -10000 100000]);
    legend('Residential','Ind & Comm','Heavy Ind','Other *');
    hold on;
    title('Forecast Total New Zealand Demand');
    ylabel('GWh');    
    drawnow;
    hold off;
end;

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