% New Zealand National and Regional 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.

% 2007 version of the model requires DemandForecastInputTablesSept2007.mat workspace to be in same directory as the script
% 2009 version requires "DF inputs.xls" to be in the 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)


% Regional allocation split between residential demand and commercial/industrial demand

% Residential demand is split on the basis of projected population for each network area (EPB) and scaled back to the national total. 
% Formula is ResDemand(FCYear) = EPBPop(FCYear)/EPBPop(BaseYear) * Nat.ResDemandPerPerson(FCYear)/Nat.ResDemandPerPerson(BaseYear) * Res.Demand(BaseYear)
% Demand is then allocated to GXP based on recent GXP demand

% Commercial/Industrial Demand is allocated to regions based on projected GDP for each region
% and scaled back to the national total
% Formula is IndCommDemand(FCYear) = GDP(FCYear)/GDP(BaseYear) * Nat.IndCommDemandperGDP(FCYear)/NatIndCommDemandperGDP(BaseYear) * IndCommDemand(BaseYear)
% Demand is then allocated to each EPB and then to GXP based on recent demand compared to total regional and network demand respectively.

% Regional allocations can be adjusted using a simple trend derived from historical regional demand (weighted relative to the population 
% and GDP based allocation using a logistics curve for the transition)
clear all
tic

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

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


% Model Switches

SingleRun = 0; % 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 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)
DisplayRegionChecks = 0; % Set this to 1 to show summary totals for regional variables
ShowGraphs = 0; % Set to 1 to show graphs
SaveRegionalMCResults = 1; % Set to 1 to save region-level Monte Carlo results to file 
SaveGXPMCResults = 0; % Set to 1 to save GXP-level Monte Carlo results to file

IncludeSpecifiedNewLoadAdjustments = 1; % set this to 1 to include new load adjustments specified in the 'NewLoad' table.
AdjustAllocByRegionalTrend = 1; % set this to 1 to adjust the regional forecast allocation using recent regional demand growth trends
ExcludeSpecHistDataFromTrend = 1; % set this to 1 to exclude specified historical load from the trend calculation (used to exclude historical Tiwai growth from the Otago/Southland trend)   
RegionalTrendAdjType = 4; % sets method of phasing out the trend adjustment 1 = reducing linear, 2 = reducing exponential, 3 = fixed, 4 = Logistic (Sigmoid). (only applies if using above regional trend adjustment) 
ExcludeNegTrendRegionFromAdj = 1; % set to 1 to exclude any regions with a negative trend from being adjusted (other than as part of scaling to NZ totals) - this is to exclude Central which has been reducing due to embedded generation growth (only applies if using above regional trend adjustment) 

UseInterRegionPopVar = 1; % set this to 1 to include uncertainty in EPB inter-area population changes 

% 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);
nReg = length(Regions);
nGXP = length(GXPnames);
nEPB = length(EPBFullNames);

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;

DirectConnectEPBStart = 50; % Position of first direct connect EPB in the EPB matrices (must all be direct connects after that point)
TiwaiEPB = 63; % Position of Tiwai direct-connect EPB in the EPB matrices

% Forecast input defined variation

FCRegionPopVarStdDev = 0.1; % Standard deviation of proportion to apply to EPB population during Monte-Carlo for modelling inter-regional population movement uncertainty
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

% Regional Allocation tables

RegionStartDate = [1998 1998 1998 1998 1998 1998 1998 1998 1998 1998 1998 1998 1998]; % for calculating trended regional values
RegionalTrendSpecPhasing = [1 .8 .6 .4 .2]; % amount of weighting applied to the trend in each year - rest of years are set to 0s 
RegionalTrendPropFactor = 0.8; % amount of reduction applied as a factor to the trend for the proportional reduction option (eg 0.8 would produce a series 1, 0.8, 0.64, 0.512 ....)
RegionalTrendFixed = 0.5; % factor applied to regional trend allocation in each year (i.e. 1 uses trend only, 0 would use POP and GDP allocation only) 
RegionTrendLogCurveDiv = 2; % Factor to apply to forecast year for calculating logistic curve value
RegionTrendLogCurveSub = 5; % Figure subtracted from divided forecast year for calculating logistic curve value
RegionTrendLogCurveDivSD = .2; % Standard deviation for logistic curve divisor for Monte Carlo
RegionTrendLogCurveSubSD = .5; % Standard deviation for logistic curve subtractor for Monte Carlo

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

% 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;
TempFCRegionGDP = FCRegionGDP;

% 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);
ForecastRegionTotalStore = zeros(ForecastPeriods,nReg,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);
IndRsquaredStore = zeros(1,MCRUNS);
CommRsquaredStore = zeros(1,MCRUNS);
IndCommRsquaredStore = zeros(1,MCRUNS);
IndCommCorrectedRsquaredStore = zeros(1,MCRUNS);
HeavyIndRsquaredStore = zeros(1,MCRUNS);
HeavyIndCorrectedRsquaredStore = zeros(1,MCRUNS);

if SaveGXPMCResults == 1
    ForecastGXPTotalStore = zeros(ForecastPeriods,nGXP,MCRUNS);
end;


% Plotting bin storage

NumberPlottingBins = 10;

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

% Precalculate recent EPB demand

EPBRecentDemand = zeros(1,nEPB);
RegionTotalDemand = zeros(1,nReg);

for gxp = 1 : nGXP;
    EPBRecentDemand(GXPEPBNum(gxp)) = EPBRecentDemand(GXPEPBNum(gxp))+ GXPRecentDemand(gxp);
end;

for epb = 2 : nEPB;
    RegionTotalDemand(EPBRegion(epb)) = RegionTotalDemand(EPBRegion(epb))+EPBRecentDemand(epb); 
end;

% Create indexes for use in regional allocations

NonDirectConnectGXPIndex = find(GXPEPBNum<DirectConnectEPBStart)';
NonTiwaiGXPIndex = find(GXPEPBNum~=TiwaiEPB)';
NonTiwaiEPBIndex = 1:nEPB; 
NonTiwaiEPBIndex(TiwaiEPB-1)=[]; 
NonTiwaiEPBIndex(1)=[]; % EPB 1 = Off grid

NonDirectConnectEPBIndex = 2:DirectConnectEPBStart-1;
 
% Start of Monte Carlo

disp(' ');
disp('Regional 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 = randn * GDPProductivityStdDev; %assumes normal distribution
        FCHouseSizeVar = randn * HouseSizeStdDev;
        FCPriceVar = randn * PriceStdDev;

        % Calculate historical model input variation for model error
        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+randn*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 forecast error
            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 Models
    % 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
    
    PropDirectConnect = 1 - sum(GXPRecentDemand(NonDirectConnectGXPIndex))./(sum(GXPRecentDemand)-GXPRecentDemand(find(GXPEPBNum==TiwaiEPB))); % Proportion direct connect load of total Domestic and Industrial/Commercial (excluding Tiwai)
    
    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];

% START REGIONAL
    
    % Set up required summary stats 
    % Scale recent demand
    
    TiwaiEPBDemand = EPBRecentDemand(TiwaiEPB);
    RegionTotalDemandlessTiwai = RegionTotalDemand - [zeros(1,length(RegionTotalDemand)-1) TiwaiEPBDemand]; % Otago/Southland in last region
    EPBRecentDemandlessTiwai = EPBRecentDemand - [zeros(1,length(EPBRecentDemand)-1) TiwaiEPBDemand]; % Tiwai last EPB
    
    ScaledEPBRecentDemandlessTiwai = EPBRecentDemandlessTiwai.*((DomDemand(length(DomDemand))+IndCommDemand(length(IndCommDemand)))./sum(EPBRecentDemandlessTiwai));
    ScaledRegionTotalDemandlessTiwai = RegionTotalDemandlessTiwai.*((DomDemand(length(DomDemand))+IndCommDemand(length(IndCommDemand)))./sum(RegionTotalDemandlessTiwai));
    
    ScaledEPBRecentDemand = ScaledEPBRecentDemandlessTiwai + [zeros(1,length(ScaledEPBRecentDemandlessTiwai)-1) TiwaiDemand(length(TiwaiDemand))];
    ScaledRegionTotalDemand = ScaledRegionTotalDemandlessTiwai + [zeros(1,length(ScaledRegionTotalDemandlessTiwai)-1) TiwaiDemand(length(TiwaiDemand))];
    % (note : scaled not matching demand because ignoring embedded vs local lines losses)
    
    % Scale regional Dom Demand and Ind Comm Demand to match EPB demand (because the breakdown is from 04/05)

    CombinedRegionIndCommDomDemand = RegionsDomDemand+RegionsIndCommDemand;
    RegionScaleFactor = ScaledRegionTotalDemand'./CombinedRegionIndCommDomDemand;
    for reg = 1:nReg;
        ScaledRegionsDomDemand = RegionsDomDemand .* RegionScaleFactor;
        ScaledRegionsIndCommDemand = RegionsIndCommDemand .* RegionScaleFactor;
    end;  
    
    % Calculate proportion of domestic load in each region (for each non direct connect EPB)
    
    RegionsIndCommDemandLessTiwai = ScaledRegionsIndCommDemand;
    RegionsIndCommDemandLessTiwai(13) = RegionsIndCommDemandLessTiwai(13) - TiwaiDemand(length(TiwaiDemand));
    
    RegionsDemandLessDC = zeros(1,nReg);
    for epb = 2: DirectConnectEPBStart-1;
        RegionsDemandLessDC(EPBRegion(epb)) = RegionsDemandLessDC(EPBRegion(epb)) + ScaledEPBRecentDemand(epb);
    end;
    
    RegionDomDemandProp = ScaledRegionsDomDemand./(RegionsDemandLessDC');
    
    % correct starting epb data to match back to total pop

    EPBPopScaleFactor = TotalNZPop(length(TotalNZPop))./sum(EPBCurrentPop(2:length(EPBCurrentPop))); % Excluding non-grid because national forecast attached to total population doing this to avoid distorting EPB forecasts
    EPBRecentPop = [0 EPBCurrentPop(2:length(EPBCurrentPop))] .*EPBPopScaleFactor;

    % scale recent region GDP figures

    RecentGDPScaleFactor = GDP(length(GDP))./sum(RecentRegionGDP);
    ScaledRecentRegionGDP = RecentRegionGDP.*RecentGDPScaleFactor;

    % Regional Forecast Allocations
    
    % Residential Allocation
    
    EPBPopGrowthProp = zeros(ForecastPeriods,nEPB);
    DomDemandPerPersonGrowthProp = zeros(1,ForecastPeriods);
    ForecastEPBDomDemand = zeros(ForecastPeriods,nEPB);
    ForecastGXPDomDemand = zeros(ForecastPeriods,nGXP);
    EPBRecentDomDemand = zeros(1,nEPB);
    
    % Scale EPB pop demand to national level pop demand forecast
    
    TotalEPBPop = sum(FCEPBPop((1:ForecastPeriods),:),2);
    OffGridPop = FCEPBPop(:,1);
    for epb = 1:length(FCEPBPop);
        ScaledFCEPBPop(:,epb)  = FCEPBPop(1:ForecastPeriods,epb).*(FCNZPOP(1:ForecastPeriods)./(TotalEPBPop-OffGridPop(1:ForecastPeriods))); % adjust for off-grid pop
    end;
    ScaledFCEPBPop(:,1)=zeros(ForecastPeriods,1); % remove the off-grid population (to keep consistent vs total demand forecasts - marginal impact - 0.1% of pop
    ScaledPopStore = ScaledFCEPBPop;
    
    % Apply inter-regional population variation and adjust region GDP for consistency 
    
    if (SingleRun ~= 1 && UseInterRegionPopVar == 1)
    
        % randomly adjust path in each region
        for epb = 1:length(EPBCurrentPop)           
            FCRegEndPopNew = (1+randn*FCRegionPopVarStdDev)* ScaledFCEPBPop(ForecastPeriods,epb);           
            FCRegPOPNew = ScaledFCEPBPop(:,epb)+(FCRegEndPopNew-ScaledFCEPBPop(ForecastPeriods,epb))*(FCYear-FCYear(1)+1)./(FCYear(ForecastPeriods)-FCYear(1)+1);
            ScaledFCEPBPop(:,epb) = FCRegPOPNew;
        end;
        % scale back to national totals again
        for year = 1:ForecastPeriods   
            ScaledFCEPBPop(year,:) = ScaledFCEPBPop(year,:) .* FCNZPOP(year)./sum(ScaledFCEPBPop(year,:),2);
        end;
        % adjust raw region GDP foreacst
        for reg = 1:nReg
            FCRegPOPnewrate = sum(ScaledFCEPBPop(2:ForecastPeriods,find(EPBRegion==reg)),2)./sum(ScaledFCEPBPop(1:ForecastPeriods-1,find(EPBRegion==reg)),2);
            FCRegPOPoldrate = sum(ScaledPopStore(2:ForecastPeriods,find(EPBRegion==reg)),2)./sum(ScaledPopStore(1:ForecastPeriods-1,find(EPBRegion==reg)),2);
            FCNZPOPVar = FCRegPOPnewrate./FCRegPOPoldrate;
            CumRegPopIncr = 1;
            for year = 2:ForecastPeriods
                CumRegPopIncr = CumRegPopIncr * FCNZPOPVar(year-1);
                FCRegionGDP(year,reg) = TempFCRegionGDP(year,reg) .* CumRegPopIncr;
            end;
        end;
    end;
        
    % Calculate factors for EPB forecasts and calc EPB level demand
      
    for epb = 2:DirectConnectEPBStart-1;
        EPBRecentDomDemand(epb) = ScaledEPBRecentDemand(epb).* RegionDomDemandProp(EPBRegion(epb));
    end;
    
    for year = 1:ForecastPeriods;
        DomDemandPerPersonGrowthProp(year) = (DomTotalForecastDemand(year)./FCNZPOP(year))./(DomDemand(EndDom)./TotalNZPop(EndDom));
        for epb = 2 : DirectConnectEPBStart-1;
            if EPBRecentPop(epb)==0; EPBPopGrowthProp(year,epb) = 0; else EPBPopGrowthProp(year,epb) = ScaledFCEPBPop(year,epb)./EPBRecentPop(epb); end;
            ForecastEPBDomDemand(year,epb) = EPBRecentDomDemand(1,epb).* EPBPopGrowthProp(year,epb).*DomDemandPerPersonGrowthProp(year);
        end;
        ForecastEPBDomDemand(year,:) = ForecastEPBDomDemand(year,:).*(DomTotalForecastDemand(year)./sum(ForecastEPBDomDemand(year,:))); % scale back to national total
    end;
    
    CheckRawEPBFCDomDemand = sum(sum(ForecastEPBDomDemand,2)-DomTotalForecastDemand);
    
    % and allocate to GXP based on proportion of recent GXP dom demand data to the relevant EPB demand
    
    for gxp = NonDirectConnectGXPIndex;
         GXPPropOfEPB = GXPRecentDemand(gxp)./EPBRecentDemand(GXPEPBNum(gxp));
         ForecastGXPDomDemand(:,gxp) = GXPPropOfEPB.*ForecastEPBDomDemand(:,GXPEPBNum(gxp));
    end;
     
    CheckGXPDomTotal=sum(sum(ForecastGXPDomDemand,2)-DomTotalForecastDemand);
     
    % and summarise back up to a regional level
    
    ForecastRegionDomDemand = zeros(ForecastPeriods,nReg);
    for reg = 1:nReg;
        ForecastRegionDomDemand(:,reg)=sum(ForecastGXPDomDemand(:,strmatch(Regions(reg),GXPregion)),2);
    end;
    CheckDomRegionTotal=sum(sum(ForecastRegionDomDemand,2)-DomTotalForecastDemand);
      
    % Commercial and Industrial Allocation
    % Allocation based on GDP/demand at a regional level, scaled to national total, and GXP demand relative to regional demand at a GXP level (excluding direct connects)

    % Regional Forecast Allocations
    
    ScaledFCRegionGDP = zeros(ForecastPeriods,nReg);
    ForecastRegionIndCommDemand = zeros(ForecastPeriods,nReg);
        
    for reg = 1:nReg;        
        ScaledFCRegionGDP(:,reg) = FCRegionGDP(1:ForecastPeriods,reg).*FCGDP./sum(FCRegionGDP(1:ForecastPeriods,:),2);    
    end;
    
    for year = 1:ForecastPeriods;
        RegionGDPGrowthProp(year,:) = ScaledFCRegionGDP(year,:)./ScaledRecentRegionGDP(:,1)';
        IndCommDemandPerGDPProp(year) = (IndCommTotalForecastDemand(year)./FCGDP(year))./(IndCommDemand(EndIndComm)./GDP(EndIndComm));
        ForecastRegionIndCommDemand(year,:) = RegionsIndCommDemandLessTiwai(:,1)'.* RegionGDPGrowthProp(year,:).*IndCommDemandPerGDPProp(year);        
        ForecastRegionIndCommDemand(year,:) = ForecastRegionIndCommDemand(year,:).*(IndCommTotalForecastDemand(year)./sum(ForecastRegionIndCommDemand(year,:))); % scale back to nat total
    end;
    
    CheckRegionIndCommTotal=sum(sum(ForecastRegionIndCommDemand,2)-IndCommTotalForecastDemand);
    
    % allocate to EPB within the region 
    
    % Recent Light Ind and Comm Demand at EPB level
    
    EPBPropOfRegion = zeros(1,nEPB);
    
    for epb = 2:nEPB;
        EPBPropOfRegion(epb)=ScaledEPBRecentDemandlessTiwai(epb)./ScaledRegionTotalDemandlessTiwai(EPBRegion(epb));
        ForecastEPBIndCommDemand(:,epb) = EPBPropOfRegion(epb).*ForecastRegionIndCommDemand(:,EPBRegion(epb));
    end;
    
    CheckEPBIndCommDemand = sum(sum(ForecastEPBIndCommDemand,2)-IndCommTotalForecastDemand);
    
    % and allocate to GXPs within the EPB (note that this subset of EPBs is larger than the residential one due to including direct connects)
    
    for gxp = NonTiwaiGXPIndex;
        GXPPropOfEPB = GXPRecentDemand(gxp)./EPBRecentDemand(GXPEPBNum(gxp));
        ForecastGXPIndCommDemand(:,gxp) = GXPPropOfEPB.*ForecastEPBIndCommDemand(:,GXPEPBNum(gxp));
    end;
    
    CheckGXPIndCommDemand = sum(sum(ForecastGXPIndCommDemand,2)-IndCommTotalForecastDemand);
        
    % Heavy Industrial Allocation
    
    % Allocate in a constant proportion to recent demand by the heavy industrials at each GXP
    
    GXPFCHeavyIndDemand = zeros(ForecastPeriods,nGXP);
    
    GXPFCHeavyIndDemand(:,find(GXPEPBNum==TiwaiEPB))= HeavyIndTotalForecastDemand;
    
    CheckGXPHeavyIndTotal=sum(sum(GXPFCHeavyIndDemand,2)-HeavyIndTotalForecastDemand);
    
    % and summarise back up to a regional level
    
    RegionalHeavyIndDemand = zeros(ForecastPeriods,nReg);
    for reg = 1:nReg;
        RegionalHeavyIndDemand(:,reg)=sum(GXPFCHeavyIndDemand(:,strmatch(Regions(reg),GXPregion))')';
    end;
    CheckHeavyIndRegionTotal=sum(sum(RegionalHeavyIndDemand,2)-HeavyIndTotalForecastDemand);
    
    % Allocate local line losses at forecast line co flat rate multiplied by the domestic and commercial/industrial load                
    
    GXPFCLinesLosses = zeros(ForecastPeriods,nGXP);
    SumOfNonDirectConnectIndCommDom = sum((ForecastGXPIndCommDemand(:,NonDirectConnectGXPIndex) + ForecastGXPDomDemand(:,NonDirectConnectGXPIndex)),2);
    for gxp = NonDirectConnectGXPIndex;
          GXPFCLinesLosses(:,gxp) =  ForecastTotalLineLosses .* (ForecastGXPIndCommDemand(:,gxp) + ForecastGXPDomDemand(:,gxp))./SumOfNonDirectConnectIndCommDom;
    end;
    CheckGXPLinesLosses = sum(sum(GXPFCLinesLosses,2)-ForecastTotalLineLosses);
      
    % Allocate embedded generation based on total demand at each grid exit point (exclude Tiwai)
    
    GXPFCEmbeddedGen = zeros(ForecastPeriods,nGXP);
  
    ScaledGXPExistingEmbeddedGWh = GXPExistingEmbeddedGWh.* ExistingEmbeddedGenerationGWh ./sum(sum(GXPExistingEmbeddedGWh')');
    for fc = 1:ForecastPeriods;
        GXPFCEmbeddedGen(fc,:)=(ForecastGXPIndCommDemand(fc,:)+ForecastGXPDomDemand(fc,:))./sum(ForecastGXPIndCommDemand(fc,:)+ForecastGXPDomDemand(fc,:)).*EmbeddedGenForecastDemand(fc);
    end;
    CheckGXPEmbeddedGen = sum(sum(GXPFCEmbeddedGen,2)-EmbeddedGenForecastDemand);
    
    % GXP Total Forecast Demand
    
    GXPTotalForecastDemand = ForecastGXPDomDemand + ForecastGXPIndCommDemand + GXPFCHeavyIndDemand + GXPFCLinesLosses - GXPFCEmbeddedGen;
    
    CheckGXPTotalDemand = sum(sum(GXPTotalForecastDemand,2)-ForecastTotalNZDemand);
    
    % Summarise GXP totals up to regional totals
    
    RegionalTotalForecastDemand = zeros(ForecastPeriods,nReg);
    for reg = 1:nReg;
        RegGXPs = strmatch(Regions(reg),GXPregion);
        RegionalTotalForecastDemand(:,reg)=sum(ForecastGXPDomDemand(:,RegGXPs),2)+sum(ForecastGXPIndCommDemand(:,RegGXPs),2)...
        +sum(GXPFCHeavyIndDemand(:,RegGXPs),2)+sum(GXPFCLinesLosses(:,RegGXPs),2)-sum(GXPFCEmbeddedGen(:,RegGXPs),2);
    end;
    CheckRegionTotal=sum(sum(RegionalTotalForecastDemand,2)-ForecastTotalNZDemand);

    % Make adjustments to the regional forecasts based on trended regional historical demand if required
    
    if AdjustAllocByRegionalTrend == 1;  
        
        RegionalTotalForecastDemandOrig = RegionalTotalForecastDemand; % store original forecasts
        
        TrendedDemandModel = zeros(2,nReg);
        FittedRegionTrendDemand = zeros(length(RegionHistYear),nReg);
        TrendEndDate = RegionHistYear(length(RegionHistYear));
        TrendStartDate  = RegionHistYear(1);
        for reg = 1:nReg;
            if ExcludeSpecHistDataFromTrend == 1; % Exclude specified historical load in the region from the trend calculation but add back the max of the specified load to the trend
                TrendedDemandModel(:,reg) = [ones(TrendEndDate-RegionStartDate(reg)+1,1) (RegionStartDate(reg):TrendEndDate)']\(RegionHistAnnualTotals(RegionStartDate(reg)-TrendStartDate+1:TrendEndDate-TrendStartDate+1, reg)-RegionExcludeFromTrend(RegionStartDate(reg)-TrendStartDate+1:TrendEndDate-TrendStartDate+1, reg));
                FittedRegionTrendDemand(RegionStartDate(reg)-TrendStartDate+1:TrendEndDate-TrendStartDate+1,reg) = [ones(TrendEndDate-RegionStartDate(reg)+1,1) (RegionStartDate(reg):TrendEndDate)']*TrendedDemandModel(:,reg) + RegionExcludeFromTrend(RegionStartDate(reg)-TrendStartDate+1:TrendEndDate-TrendStartDate+1, reg);
            else                        
                TrendedDemandModel(:,reg) = [ones(TrendEndDate-RegionStartDate(reg)+1,1) (RegionStartDate(reg):TrendEndDate)']\RegionHistAnnualTotals(RegionStartDate(reg)-TrendStartDate+1:TrendEndDate-TrendStartDate+1, reg);
                FittedRegionTrendDemand(RegionStartDate(reg)-TrendStartDate+1:TrendEndDate-TrendStartDate+1,reg) = [ones(TrendEndDate-RegionStartDate(reg)+1,1) (RegionStartDate(reg):TrendEndDate)']*TrendedDemandModel(:,reg);
            end;   
        end;
        
        RawRegionalTrendForecast = zeros(ForecastPeriods,nReg);
        for reg = 1:nReg;
            RawRegionalTrendForecast(:,reg) = [ones(ForecastPeriods,1) FCYear]*TrendedDemandModel(:,reg) + ((ExcludeSpecHistDataFromTrend == 1).* + max(RegionExcludeFromTrend(:,reg)));
        end;

        if ExcludeNegTrendRegionFromAdj == 1
            TrendRegions = TrendedDemandModel(2,:)>0; % identifies regions with a positive trend slope
        else
            TrendRegions = ones(1,nReg);
        end;
        
        % multiply allocated forecasts by calculated ratio
        
        TrendWeighting = zeros(ForecastPeriods,1);
        for AdjYear = 1:ForecastPeriods            
            switch RegionalTrendAdjType
                case 1 % reducing linear or specified
                    if AdjYear <= length(RegionalTrendSpecPhasing)
                        TrendWeighting(AdjYear)=RegionalTrendSpecPhasing(AdjYear);
                    else
                        TrendWeighting(AdjYear)=0;
                    end;
                case 2 % reducing exponential                    
                    TrendWeighting(AdjYear) = RegionalTrendPropFactor^(AdjYear-1);                        
                case 3 % fixed 
                    TrendWeighting(AdjYear) = RegionalTrendFixed;
                case 4 % logistic curve
                    if SingleRun ~= 1
                        TrendDivVar = randn * RegionTrendLogCurveDivSD;
                        TrendSubVar = randn * RegionTrendLogCurveSubSD;
                    else
                        TrendDivVar = 0;
                        TrendSubVar = 0;
                    end;              
                    TrendWeighting(AdjYear) = 1-(1/(1+exp(-(AdjYear/(RegionTrendLogCurveDiv+TrendDivVar)-(RegionTrendLogCurveSub+TrendSubVar)))));                    
                otherwise
                    error('Trend adjustment type incorrect');
            end;
            RegionalTotalForecastDemand(AdjYear,TrendRegions) = (RegionalTotalForecastDemand(AdjYear,TrendRegions).* (1-TrendWeighting(AdjYear))) + (RawRegionalTrendForecast(AdjYear,TrendRegions).* TrendWeighting(AdjYear));
        end;
             
        % scale back to national totals again
        
        for AdjYear = 1:ForecastPeriods
            RegionalTotalForecastDemand(AdjYear,:) = RegionalTotalForecastDemand(AdjYear,:).* ForecastTotalNZDemand(AdjYear)./sum(RegionalTotalForecastDemand(AdjYear,:));
        end;
                
        % scale GXP totals to match regions

        for reg = 1:nReg        
            RegGXPs = strmatch(Regions(reg),GXPregion);
            for RegGXP = 1:length(RegGXPs)
                AdjGXP = RegGXPs(RegGXP);
                ForecastGXPDomDemand(:,AdjGXP) = ForecastGXPDomDemand(:,AdjGXP).* RegionalTotalForecastDemand(:,reg)./ RegionalTotalForecastDemandOrig(:,reg);
                ForecastGXPIndCommDemand(:,AdjGXP) = ForecastGXPIndCommDemand(:,AdjGXP).* RegionalTotalForecastDemand(:,reg)./ RegionalTotalForecastDemandOrig(:,reg);
                GXPFCHeavyIndDemand(:,AdjGXP) = GXPFCHeavyIndDemand(:,AdjGXP) .* RegionalTotalForecastDemand(:,reg)./ RegionalTotalForecastDemandOrig(:,reg);
                GXPFCLinesLosses(:,AdjGXP) = GXPFCLinesLosses(:,AdjGXP) .* RegionalTotalForecastDemand(:,reg)./ RegionalTotalForecastDemandOrig(:,reg);
                GXPFCEmbeddedGen(:,AdjGXP) = GXPFCEmbeddedGen(:,AdjGXP) .* RegionalTotalForecastDemand(:,reg)./ RegionalTotalForecastDemandOrig(:,reg);
            end;
        end;      
    
        GXPTotalForecastDemand = ForecastGXPDomDemand + ForecastGXPIndCommDemand + GXPFCHeavyIndDemand + GXPFCLinesLosses - GXPFCEmbeddedGen;
        
        CheckAdjGXPs = sum(sum(ForecastGXPDomDemand(:,:))+sum(ForecastGXPIndCommDemand(:,:))+sum(GXPFCHeavyIndDemand(:,:))+sum(GXPFCLinesLosses(:,:))-sum(GXPFCEmbeddedGen(:,:)))-sum(ForecastTotalNZDemand);
        CheckAdjGXPTotal = sum(sum(GXPTotalForecastDemand,2)-ForecastTotalNZDemand);
        
    end;
    
    % Make any specified adjustments to forecasts
    % Known major new loads (anything over 5% of existing GXP load as a guideline) should be included in the "NewLoad" table 
    % in the Matlab workspace - the fields in the table are: NewLoadName, NewLoadBusName, NewLoadMW, NewLoadGWh, NewLoadStartYear, NewLoadEndYear, 
    % MWLoadAdditiveFlag GWhLoadAdditiveFlag GXPLoadDiversity, ExtremeSummerRegionLoadDiversity, SummerRegionLoadDiversity,WinterRegionLoadDiversity, 
    % ExtremeSummerIslandLoadDiversity, SummerIslandLoadDiversity , WinterIslandLoadDiversity, ExtremeSummerIslandLightLoadDiversity, SummerIslandLightLoadDiversity, WinterIslandLightLoadDiversity
    % Where the new load has no scheduled 'end date' use 2099 (any date past the last forecast year will have the same effect)
    % Generally load should be non-additive to the total national demand forecast - i.e. adjustments are proportionately spread and 
    % subtracted off all other GXPs so that the national total remains the same but the NewLoadAdditiveFlag can be set to 1 to force the load 
    % to be added to the total national forecast if desired.
    % New Load MW and the diversity figures are used for peak forecasting. Diversity is assumed at 100% by default - ie in the absense of figures it is 
    % assumed that 100% of the MW load will impact on peak GXP, regional, and island load. 
    % Avoidance of double counting of growth through both GWh and MW changes is handled in the ADMD script  

    GXPFCNewLoadAdj = zeros(ForecastPeriods,nGXP);
        
    if IncludeSpecifiedNewLoadAdjustments == 1; % only do the new load calculations if the switch is 'On'
    
        % calculate total GXP demand less Heavy Industrial

        GXPTotalDemandLessHeavyInd = ForecastGXPDomDemand + ForecastGXPIndCommDemand + GXPFCLinesLosses - GXPFCEmbeddedGen;

        % accumulate specified new loads

        if exist('NewLoad')==1 && not(isempty(NewLoadText{1,1}));
            NewLoadDim = size(NewLoad);
            for NewLoadNum = 1:NewLoadDim(1);
                NewLoadBusName = NewLoad{NewLoadNum,2};
                NewLoadMWh = NewLoad{NewLoadNum,4};
                NewLoadStartYear = NewLoad{NewLoadNum,5};
                NewLoadEndYear = NewLoad{NewLoadNum,6};
                NewLoadAdditiveFlag = NewLoad{NewLoadNum,8};
                LastActualYear = Year(length(Year));
                AdjGXP = strmatch(NewLoadBusName,GXPnames);
                if not(isempty(AdjGXP))
                    if NewLoadStartYear<LastActualYear+1;
                        NewLoadStartYear = LastActualYear+1;
                        disp(['Warning - new load start date is earlier than the first forecast period for ' NewLoad{NewLoadNum,1}]);
                    end;
                    if NewLoadEndYear>LastActualYear+ForecastPeriods;
                        NewLoadEndYear=LastActualYear+ForecastPeriods;
                    end;
                    if NewLoadEndYear>=NewLoadStartYear;
                        for AdjYear = 1:ForecastPeriods;
                            if (AdjYear + LastActualYear >= NewLoadStartYear ) & (AdjYear + LastActualYear <= NewLoadEndYear)
                                % Add load to the specified bus
                                GXPFCNewLoadAdj(AdjYear,AdjGXP) = GXPFCNewLoadAdj(AdjYear,AdjGXP) + NewLoadMWh;
                                % If not additive, spread and subtract the load from all other buses
                                if NewLoadAdditiveFlag ~= 1;
                                    GXPAdjProportion = NewLoadMWh./(sum(GXPTotalDemandLessHeavyInd(AdjYear,:))-GXPTotalDemandLessHeavyInd(AdjYear,AdjGXP));
                                    for gxp = 1:nGXP;
                                        if gxp ~= AdjGXP
                                            GXPFCNewLoadAdj(AdjYear,gxp) = GXPFCNewLoadAdj(AdjYear,gxp) - (GXPTotalDemandLessHeavyInd(AdjYear,gxp).* GXPAdjProportion);
                                        end;
                                    end;
                                end;
                            end;
                        end;
                    else
                        disp(['Warning - new load start date is later than the new load end date for ' NewLoad{NewLoadNum,1}]);
                    end;
                else
                    disp(['Warning - new load bus name does not match existing bus names for ' NewLoad{NewLoadNum,1}]);
                end;
            end;
        end;    
    end;
    
    % adjust region totals
    
    for reg = 1:nReg;
        RegionalTotalForecastDemand(:,reg)=RegionalTotalForecastDemand(:,reg)+sum(GXPFCNewLoadAdj(:,strmatch(Regions(reg),GXPregion)),2);
    end;
           
    GXPTotalForecastDemandAdj = GXPTotalForecastDemand + GXPFCNewLoadAdj;
    
    CheckNewLoadAdj = sum(sum(RegionalTotalForecastDemand')'-ForecastTotalNZDemand);
    CheckGXPTotalDemandAdj = sum(sum(GXPTotalForecastDemandAdj,2)-ForecastTotalNZDemand);
    
% END REGIONAL
    
    % 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;
    ForecastRegionTotalStore(:,:,MCCycle) = RegionalTotalForecastDemand;
    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;
    
    if SaveGXPMCResults == 1
        ForecastGXPTotalStore(:,:,MCCycle) = GXPTotalForecastDemandAdj;
    end;
    
    for i=1:ForecastPeriods;
        TotalDemandBin(i,fix(ForecastTotalNZDemand(i)/(120000/NumberPlottingBins)))=TotalDemandBin(i,fix(ForecastTotalNZDemand(i)/(120000/NumberPlottingBins)))+1;
    end

end % of MonteCarlo

%xlswrite code
ExcelOutput1 = [FCYear,DomTotalForecastDemand, IndCommTotalForecastDemand,HeavyIndTotalForecastDemand,ForecastTotalLineLosses,EmbeddedGenForecastDemand,ForecastTotalNZDemand];
ExcelOutput2 = [RegionalTotalForecastDemand];
ExcelOutput3 = [GXPTotalForecastDemandAdj];

xlswrite(strcat(OutputDir,OutputFN),ExcelOutput1,'RegForecast', 'a10');
xlswrite(strcat(OutputDir,OutputFN),ExcelOutput2,'RegForecast', 'i10');
xlswrite(strcat(OutputDir,OutputFN),ExcelOutput3,'RegForecast', 'w10');


% RegionCheck - for checking individual regions (test purposes)
if DisplayRegionChecks == 1;
    for reg  = 1:nReg
        disp(cellstr(Regions(reg)));
        disp('RegHistTotal   RegRecentTotal   ScaledHistReg   FirstDomFC');
        disp('EstRegionDom   FirstIndCommFC   EstRegionInd    FirstLineLoss   FirstEMB');
        disp([num2str(RegionHistAnnualTotals(size(RegionHistAnnualTotals,1),reg)), '           ', num2str(RegionTotalDemand(reg)), '         ', num2str(ScaledRegionTotalDemandlessTiwai(1,reg)), '       ', num2str(sum(ForecastGXPDomDemand(1,strmatch(Regions(reg),GXPregion)),2))]);
        disp([num2str(ScaledRegionsDomDemand(reg)),'      ', num2str(sum(ForecastGXPIndCommDemand(1,strmatch(Regions(reg),GXPregion)),2)), '         ', num2str(ScaledRegionsIndCommDemand(reg)),  '        ',  num2str(sum(GXPFCLinesLosses(1,strmatch(Regions(reg),GXPregion)),2)), '        ', num2str(sum(GXPFCEmbeddedGen(1,strmatch(Regions(reg),GXPregion)),2))]);
        disp('          ');
    end;
end;

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);
IndCommForecastLowerPercentile = zeros(1,ForecastPeriods);
IndCommForecastUpperPercentile = zeros(1,ForecastPeriods);
IndCommForecastMedian = 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 = 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));

    for reg = 1:nReg;
        SortedForecast = zeros(ForecastPeriods,MCRUNS);
        for i = 1:ForecastPeriods;
            SortedForecast(i,:)=(ForecastRegionTotalStore(i,reg,:));
        end;
        for i=1:ForecastPeriods;
            SortedForecast(i,:)=sort(SortedForecast(i,:));
        end;
        RegionTotalForecastLowerPercentile(reg,:) = SortedForecast(:,fix(ForecastPercentile.*MCRUNS));
        RegionTotalForecastUpperPercentile(reg,:)  = SortedForecast(:,MCRUNS-fix(ForecastPercentile.*MCRUNS));
        RegionTotalForecastMedian(reg,:)  = SortedForecast(:,fix(MCRUNS./2));
    end;
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');
    disp(' Figure 23 : GXP level residential demand');
    disp(' Figure 24 : GXP level industrial and commercial demand');
    disp(' Figure 25 : GXP level heavy industrial demand');
    disp(' Figure 26 : Regional residential demand');
    disp(' Figure 27 : Regional industrial and commercial demand');
    disp(' Figure 28 : Regional heavy industrial demand');
    disp(' Figure 29 : Regional total demand');
    disp(' Figure 30 : Total demand Upper North Island');
    disp(' Figure 31 : Total demand Upper South Island');
    disp(' Figure 32 : Total demand North Isthmus');
    disp(' Figure 33 : Total demand Auckland');
    disp(' Figure 34 : Total demand Waikato');
    disp(' Figure 35 : Total demand Bay of Plenty');
    disp(' Figure 36 : Total demand Hawkes Bay');
    disp(' Figure 37 : Total demand Central');
    disp(' Figure 38 : Total demand Taranaki');
    disp(' Figure 39 : Total demand Wellington');
    disp(' Figure 40 : Total demand Nelson / Marlborough');
    disp(' Figure 41 : Total demand West Coast');
    disp(' Figure 42 : Total demand Canterbury');
    disp(' Figure 43 : Total demand South Canterbury');
    disp(' Figure 44 : Total demand Otago / Southland');
    
    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 Residential Demand');
    ylabel('GWh');
    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');
    ylabel('GWh');
    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');
    ylabel('GWh');
    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');
    ylabel('GWh');
    drawnow;
    hold off;
    
    figure(8);
    hist(DomRsquaredStore);
    title('Residential 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 Residential 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;
    ylabel('GWh');
    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;
    ylabel('GWh');
    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;
    ylabel('GWh');
    drawnow;
    hold off;

    % The graph produced in figure 14 has embedded generation as a negative so that the total demand is correct. 
    % The embedded generation 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');
    ylabel('GWh');
    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;

    figure(23);
    ForecastGXPDomGraph  = surf(1:nGXP, FCYear, ForecastGXPDomDemand);
    title('GXP Residential Forecast Demand');
    drawnow;

    figure(24);
    ForecastGXPIndCommGraph = surf(1:nGXP, FCYear, ForecastGXPIndCommDemand);
    title('GXP Commercial and Industrial Forecast Demand');
    drawnow;

    figure(25);
    GXPHeavyIndForecastGraph = surf(1:nGXP, FCYear, GXPFCHeavyIndDemand);
    title('GXP Heavy Industrial Forecast Demand');
    drawnow;

    figure(26);
    ForecastRegionalDomGraph = surf(1:nReg, FCYear, ForecastRegionDomDemand);
    title('Regional Residential Forecast Demand');
    drawnow;

    figure(27);
    ForecastRegionalIndCommGraph = surf(1:nReg, FCYear, ForecastRegionIndCommDemand);
    title('Regional Commercial and Industrial Forecast Demand');
    drawnow;

    figure(28);
    RegionalHeavyIndForecastGraph = surf(1:nReg, FCYear, RegionalHeavyIndDemand);
    title('Regional Heavy Industrial Forecast Demand');
    drawnow;

    figure(29);
    RegionalTotalDemandGraph = surf(1:nReg, FCYear, RegionalTotalForecastDemand);
    axis([1 nReg FCYear(1) FCYear(ForecastPeriods) 0 20000 0 15000]);
    title('Regional Total Forecast Demand');
    drawnow;

    figure(30);
    % Plot of Combined Auckland and North Isthmus region
    NorthRegionPlotData=zeros(ForecastPeriods,MCRUNS);
    NorthRegionPlotData(:,:)=ForecastRegionTotalStore(:,1,:)+ForecastRegionTotalStore(:,2,:);
    % and plot
    RegionalResultSpread=plot(FCYear(1:ForecastPeriods),NorthRegionPlotData);
    hold on;
    axis([1998 FCYear(ForecastPeriods) 0 45000]);
    title('Auckland and North Isthmus Forecast Growth');
    HistoricalAucklandNorthIsthmusGWH = RegionHistAnnualTotals(:,1)+RegionHistAnnualTotals(:,2);
    HistAuckPlot = plot(RegionHistYear,HistoricalAucklandNorthIsthmusGWH, 'k-');
    set(HistAuckPlot,'LineWidth',3);
    % Calc percentiles for combined region 
    if SingleRun~=1;
        SortedForecast=NorthRegionPlotData;
        for i=1:ForecastPeriods;
            SortedForecast(i,:)=sort(SortedForecast(i,:));
        end;
        NorthRegionLowerPercentile = SortedForecast(:,fix(ForecastPercentile.*MCRUNS));
        NorthRegionUpperPercentile = SortedForecast(:,MCRUNS-fix(ForecastPercentile.*MCRUNS));
        NorthRegionMedian  = SortedForecast(:,fix(MCRUNS./2));
        UpperRegionPCGraph = plot(FCYear(1:ForecastPeriods),NorthRegionUpperPercentile,'k-');
        LowerRegionPCGraph = plot(FCYear(1:ForecastPeriods),NorthRegionLowerPercentile,'k-');
        MedianRegionForecastGraph = plot(FCYear(1:ForecastPeriods),NorthRegionMedian,'k-');
        set(UpperRegionPCGraph,'LineWidth',3);
        set(LowerRegionPCGraph,'LineWidth',3);
        set(MedianRegionForecastGraph,'LineWidth',3);
    else
        NorthRegionMedian = NorthRegionPlotData(1,:);
    end;
    drawnow;
    hold off;

    figure(31);
    % Plot of Combined Northern-South Island Region
    NorthernSouthPlotData=zeros(ForecastPeriods,MCRUNS);
    NorthernSouthPlotData(:,:)=ForecastRegionTotalStore(:,9,:)+ForecastRegionTotalStore(:,10,:)+ForecastRegionTotalStore(:,11,:);    
    % and plot
    RegionalResultSpread=plot(FCYear(1:ForecastPeriods),NorthernSouthPlotData);
    hold on;
    axis([1998 FCYear(ForecastPeriods) 0 20000]);
    title('Canterbury, Nelson-Marlborough and West Coast Forecast Growth');
    HistoricalNorthernSouthIslandGWH = RegionHistAnnualTotals(:,9)+RegionHistAnnualTotals(:,10)+RegionHistAnnualTotals(:,11);
    HistNSPlot = plot(RegionHistYear,HistoricalNorthernSouthIslandGWH, 'k-');
    set(HistNSPlot,'LineWidth',3);
    if SingleRun~=1;
        % Calc percentiles for combined region 
        SortedForecast=NorthernSouthPlotData;
        for i=1:ForecastPeriods;
            SortedForecast(i,:)=sort(SortedForecast(i,:));
        end;
        NorthernSouthLowerPercentile = SortedForecast(:,fix(ForecastPercentile.*MCRUNS));
        NorthernSouthUpperPercentile = SortedForecast(:,MCRUNS-fix(ForecastPercentile.*MCRUNS));
        NorthernSouthMedian  = SortedForecast(:,fix(MCRUNS./2));
        NorthernSouthUpperPCGraph = plot(FCYear(1:ForecastPeriods),NorthernSouthUpperPercentile,'k-');
        NorthernSouthLowerPCGraph = plot(FCYear(1:ForecastPeriods),NorthernSouthLowerPercentile,'k-');
        NorthernSouthMedianGraph = plot(FCYear(1:ForecastPeriods),NorthernSouthMedian,'k-');
        set(NorthernSouthUpperPCGraph,'LineWidth',3);
        set(NorthernSouthLowerPCGraph,'LineWidth',3);
        set(NorthernSouthMedianGraph,'LineWidth',3);
    else
        NorthernSouthMedian = NorthernSouthPlotData(1,:);
    end;
    drawnow;
    hold off;

    figure(32);
    % Plot of North Isthmus
    NorthIsthmusPlotData=zeros(ForecastPeriods,MCRUNS);
    NorthIsthmusPlotData(:,:)=ForecastRegionTotalStore(:,1,:);  
    HistoricalNorthIsthmusGWH = RegionHistAnnualTotals(:,1);
    % and plot
    RegionalResultSpread=plot(FCYear(1:ForecastPeriods),NorthIsthmusPlotData);
    hold on;
    axis([1998 FCYear(ForecastPeriods) 0 15000]);
    title('NorthIsthmus Forecast Growth');
    HistNorthIsthmusPlot = plot(RegionHistYear,HistoricalNorthIsthmusGWH, 'k-');
    set(HistNorthIsthmusPlot,'LineWidth',3);
    if SingleRun~=1;
        NorthIsthmusUpperPCGraph = plot(FCYear(1:ForecastPeriods),RegionTotalForecastUpperPercentile(1,:),'k-');
        NorthIsthmusLowerPCGraph = plot(FCYear(1:ForecastPeriods),RegionTotalForecastLowerPercentile(1,:),'k-');
        NorthIsthmusMedianGraph = plot(FCYear(1:ForecastPeriods),RegionTotalForecastMedian(1,:),'k-');
        set(NorthIsthmusUpperPCGraph,'LineWidth',3);
        set(NorthIsthmusLowerPCGraph,'LineWidth',3);
        set(NorthIsthmusMedianGraph,'LineWidth',3);
    else
        NorthIsthmusMedian = NorthIsthmusPlotData(1,:);
    end;
    drawnow;
    hold off;
 
    figure(33);
    % Plot of Auckland
    AKPlotData=zeros(ForecastPeriods,MCRUNS);
    AKPlotData(:,:)=ForecastRegionTotalStore(:,2,:);
    HistoricalAKGWH = RegionHistAnnualTotals(:,2);
    % and plot
    RegionalResultSpread=plot(FCYear(1:ForecastPeriods),AKPlotData);
    hold on;
    axis([1998 FCYear(ForecastPeriods) 0 30000]);
    title('Auckland Forecast Growth');
    HistAKPlot = plot(RegionHistYear,HistoricalAKGWH, 'k-');
    set(HistAKPlot,'LineWidth',3);
    if SingleRun~=1;
        AKUpperPCGraph = plot(FCYear(1:ForecastPeriods),RegionTotalForecastUpperPercentile(2,:),'k-');
        AKLowerPCGraph = plot(FCYear(1:ForecastPeriods),RegionTotalForecastLowerPercentile(2,:),'k-');
        AKMedianGraph = plot(FCYear(1:ForecastPeriods),RegionTotalForecastMedian(2,:),'k-');
        set(AKUpperPCGraph,'LineWidth',3);
        set(AKLowerPCGraph,'LineWidth',3);
        set(AKMedianGraph,'LineWidth',3);
    else
        AKMedian = AKPlotData(1,:);
    end;
    drawnow;
    hold off;

    figure(34);
    % Plot of Waikato
    WaikatoPlotData=zeros(ForecastPeriods,MCRUNS);
    WaikatoPlotData(:,:)=ForecastRegionTotalStore(:,3,:);  
    HistoricalWaikatoGWH = RegionHistAnnualTotals(:,3);
    % and plot
    RegionalResultSpread=plot(FCYear(1:ForecastPeriods),WaikatoPlotData);
    hold on;
    axis([1998 FCYear(ForecastPeriods) 0 10000]);
    title('Waikato Forecast Growth');
    HistWaikatoPlot = plot(RegionHistYear,HistoricalWaikatoGWH, 'k-');
    set(HistWaikatoPlot,'LineWidth',3);
    if SingleRun~=1;
        WaikatoUpperPCGraph = plot(FCYear(1:ForecastPeriods),RegionTotalForecastUpperPercentile(3,:),'k-');
        WaikatoLowerPCGraph = plot(FCYear(1:ForecastPeriods),RegionTotalForecastLowerPercentile(3,:),'k-');
        WaikatoMedianGraph = plot(FCYear(1:ForecastPeriods),RegionTotalForecastMedian(3,:),'k-');
        set(WaikatoUpperPCGraph,'LineWidth',3);
        set(WaikatoLowerPCGraph,'LineWidth',3);
        set(WaikatoMedianGraph,'LineWidth',3);
    else
        WaikatoMedian = WaikatoPlotData(1,:);
    end;
    drawnow;
    hold off;

    figure(35);
    % Plot of BOP
    BOPPlotData=zeros(ForecastPeriods,MCRUNS);
    BOPPlotData(:,:)=ForecastRegionTotalStore(:,4,:);
    HistoricalBOPGWH = RegionHistAnnualTotals(:,4);
    % and plot
    RegionalResultSpread=plot(FCYear(1:ForecastPeriods),BOPPlotData);
    hold on;
    axis([1998 FCYear(ForecastPeriods) 0 10000]);
    title('BOP Forecast Growth');
    HistBOPPlot = plot(RegionHistYear,HistoricalBOPGWH, 'k-');
    set(HistBOPPlot,'LineWidth',3);
    if SingleRun~=1;
        BOPUpperPCGraph = plot(FCYear(1:ForecastPeriods),RegionTotalForecastUpperPercentile(4,:),'k-');
        BOPLowerPCGraph = plot(FCYear(1:ForecastPeriods),RegionTotalForecastLowerPercentile(4,:),'k-');
        BOPMedianGraph = plot(FCYear(1:ForecastPeriods),RegionTotalForecastMedian(4,:),'k-');
        set(BOPUpperPCGraph,'LineWidth',3);
        set(BOPLowerPCGraph,'LineWidth',3);
        set(BOPMedianGraph,'LineWidth',3);
    else
        BOPMedian = BOPPlotData(1,:);
    end;
    drawnow;
    hold off;

    
    figure(36);
    % Plot of HawkesBay
    HawkesBayPlotData=zeros(ForecastPeriods,MCRUNS);
    HawkesBayPlotData(:,:)=ForecastRegionTotalStore(:,5,:);  
    HistoricalHawkesBayGWH = RegionHistAnnualTotals(:,5);
    % and plot
    RegionalResultSpread=plot(FCYear(1:ForecastPeriods),HawkesBayPlotData);
    hold on;
    axis([1998 FCYear(ForecastPeriods) 0 5000]);
    title('HawkesBay Forecast Growth');
    HistHawkesBayPlot = plot(RegionHistYear,HistoricalHawkesBayGWH, 'k-');
    set(HistHawkesBayPlot,'LineWidth',3);
    if SingleRun~=1;
        HawkesBayUpperPCGraph = plot(FCYear(1:ForecastPeriods),RegionTotalForecastUpperPercentile(5,:),'k-');
        HawkesBayLowerPCGraph = plot(FCYear(1:ForecastPeriods),RegionTotalForecastLowerPercentile(5,:),'k-');
        HawkesBayMedianGraph = plot(FCYear(1:ForecastPeriods),RegionTotalForecastMedian(5,:),'k-');
        set(HawkesBayUpperPCGraph,'LineWidth',3);
        set(HawkesBayLowerPCGraph,'LineWidth',3);
        set(HawkesBayMedianGraph,'LineWidth',3);
    else
        HawkesBayMedian = HawkesBayPlotData(1,:);
    end;
    drawnow;
    hold off;

    figure(37);
    % Plot of Central
    CentralPlotData=zeros(ForecastPeriods,MCRUNS);
    CentralPlotData(:,:)=ForecastRegionTotalStore(:,6,:);  
    HistoricalCentralGWH = RegionHistAnnualTotals(:,6);
    % and plot
    RegionalResultSpread=plot(FCYear(1:ForecastPeriods),CentralPlotData);
    hold on;
    axis([1998 FCYear(ForecastPeriods) 0 5000]);
    title('Central Forecast Growth');
    HistCentralPlot = plot(RegionHistYear,HistoricalCentralGWH, 'k-');
    set(HistCentralPlot,'LineWidth',3);
    if SingleRun~=1;
        CentralUpperPCGraph = plot(FCYear(1:ForecastPeriods),RegionTotalForecastUpperPercentile(6,:),'k-');
        CentralLowerPCGraph = plot(FCYear(1:ForecastPeriods),RegionTotalForecastLowerPercentile(6,:),'k-');
        CentralMedianGraph = plot(FCYear(1:ForecastPeriods),RegionTotalForecastMedian(6,:),'k-');
        set(CentralUpperPCGraph,'LineWidth',3);
        set(CentralLowerPCGraph,'LineWidth',3);
        set(CentralMedianGraph,'LineWidth',3);
    else
        CentralMedian = CentralPlotData(1,:);
    end;
    drawnow;
    hold off;    
    
    figure(38);
    % Plot of Taranaki
    TaranakiPlotData=zeros(ForecastPeriods,MCRUNS);
    TaranakiPlotData(:,:)=ForecastRegionTotalStore(:,7,:);
    HistoricalTaranakiGWH = RegionHistAnnualTotals(:,7);
    % and plot
    RegionalResultSpread=plot(FCYear(1:ForecastPeriods),TaranakiPlotData);
    hold on;
    axis([1998 FCYear(ForecastPeriods) 0 2500]);
    title('Taranaki Forecast Growth');
    HistTaranakiPlot = plot(RegionHistYear,HistoricalTaranakiGWH, 'k-');
    set(HistTaranakiPlot,'LineWidth',3);
    if SingleRun~=1;
        TaranakiUpperPCGraph = plot(FCYear(1:ForecastPeriods),RegionTotalForecastUpperPercentile(7,:),'k-');
        TaranakiLowerPCGraph = plot(FCYear(1:ForecastPeriods),RegionTotalForecastLowerPercentile(7,:),'k-');
        TaranakiMedianGraph = plot(FCYear(1:ForecastPeriods),RegionTotalForecastMedian(7,:),'k-');
        set(TaranakiUpperPCGraph,'LineWidth',3);
        set(TaranakiLowerPCGraph,'LineWidth',3);
        set(TaranakiMedianGraph,'LineWidth',3);
    else
        TaranakiMedian = TaranakiPlotData(1,:);
    end;
    drawnow;
    hold off;

    figure(39);
    % Plot of Wellington
    WellingtonPlotData=zeros(ForecastPeriods,MCRUNS);
    WellingtonPlotData(:,:)=ForecastRegionTotalStore(:,8,:);  
    HistoricalWellingtonGWH = RegionHistAnnualTotals(:,8);
    % and plot
    RegionalResultSpread=plot(FCYear(1:ForecastPeriods),WellingtonPlotData);
    hold on;
    axis([1998 FCYear(ForecastPeriods) 0 10000]);
    title('Wellington Forecast Growth');
    HistWellingtonPlot = plot(RegionHistYear,HistoricalWellingtonGWH, 'k-');
    set(HistWellingtonPlot,'LineWidth',3);
    if SingleRun~=1;        
        WellingtonUpperPCGraph = plot(FCYear(1:ForecastPeriods),RegionTotalForecastUpperPercentile(8,:),'k-');
        WellingtonLowerPCGraph = plot(FCYear(1:ForecastPeriods),RegionTotalForecastLowerPercentile(8,:),'k-');
        WellingtonMedianGraph = plot(FCYear(1:ForecastPeriods),RegionTotalForecastMedian(8,:),'k-');
        set(WellingtonUpperPCGraph,'LineWidth',3);
        set(WellingtonLowerPCGraph,'LineWidth',3);
        set(WellingtonMedianGraph,'LineWidth',3);
    else
        WellingtonMedian = WellingtonPlotData(1,:);
    end;
    drawnow;
    hold off;
    
    figure(40);
    % Plot of Nelson Marlborough
    NelsonMarlPlotData=zeros(ForecastPeriods,MCRUNS);
    NelsonMarlPlotData(:,:)=ForecastRegionTotalStore(:,9,:);  
    HistoricalNelsonMarlGWH = RegionHistAnnualTotals(:,9);
    % and plot
    RegionalResultSpread=plot(FCYear(1:ForecastPeriods),NelsonMarlPlotData);
    hold on;
    axis([1998 FCYear(ForecastPeriods) 0 4000]);
    title('Nelson Marlborough Forecast Growth');
    HistNelsonMarlPlot = plot(RegionHistYear,HistoricalNelsonMarlGWH, 'k-');
    set(HistNelsonMarlPlot,'LineWidth',3);
    if SingleRun~=1;
        NelsonMarlUpperPCGraph = plot(FCYear(1:ForecastPeriods),RegionTotalForecastUpperPercentile(9,:),'k-');
        NelsonMarlLowerPCGraph = plot(FCYear(1:ForecastPeriods),RegionTotalForecastLowerPercentile(9,:),'k-');
        NelsonMarlMedianGraph = plot(FCYear(1:ForecastPeriods),RegionTotalForecastMedian(9,:),'k-');
        set(NelsonMarlUpperPCGraph,'LineWidth',3);
        set(NelsonMarlLowerPCGraph,'LineWidth',3);
        set(NelsonMarlMedianGraph,'LineWidth',3);
    else
        NelsonMarlMedian = NelsonMarlPlotData(1,:);
    end;
    drawnow;
    hold off;
    
    figure(41);
    % Plot of WestCoast
    WestCoastPlotData=zeros(ForecastPeriods,MCRUNS);
    WestCoastPlotData(:,:)=ForecastRegionTotalStore(:,10,:);
    HistoricalWestCoastGWH = RegionHistAnnualTotals(:,10);
    % and plot
    RegionalResultSpread=plot(FCYear(1:ForecastPeriods),WestCoastPlotData);
    hold on;
    axis([1998 FCYear(ForecastPeriods) 0 1000]);
    title('WestCoast Forecast Growth');
    HistWestCoastPlot = plot(RegionHistYear,HistoricalWestCoastGWH, 'k-');
    set(HistWestCoastPlot,'LineWidth',3);
    if SingleRun~=1;
        WestCoastUpperPCGraph = plot(FCYear(1:ForecastPeriods),RegionTotalForecastUpperPercentile(10,:),'k-');
        WestCoastLowerPCGraph = plot(FCYear(1:ForecastPeriods),RegionTotalForecastLowerPercentile(10,:),'k-');
        WestCoastMedianGraph = plot(FCYear(1:ForecastPeriods),RegionTotalForecastMedian(10,:),'k-');
        set(WestCoastUpperPCGraph,'LineWidth',3);
        set(WestCoastLowerPCGraph,'LineWidth',3);
        set(WestCoastMedianGraph,'LineWidth',3);
    else
        WestCoastMedian = WestCoastPlotData(1,:);
    end;
    drawnow;
    hold off;
    
    figure(42);
    % Plot of Canterbury
    CanterburyPlotData=zeros(ForecastPeriods,MCRUNS);
    CanterburyPlotData(:,:)=ForecastRegionTotalStore(:,11,:);  
    HistoricalCanterburyGWH = RegionHistAnnualTotals(:,11);
    % and plot
    RegionalResultSpread=plot(FCYear(1:ForecastPeriods),CanterburyPlotData);
    hold on;
    axis([1998 FCYear(ForecastPeriods) 0 15000]);
    title('Canterbury Forecast Growth');
    HistCanterburyPlot = plot(RegionHistYear,HistoricalCanterburyGWH, 'k-');
    set(HistCanterburyPlot,'LineWidth',3);
    if SingleRun~=1;
        CanterburyUpperPCGraph = plot(FCYear(1:ForecastPeriods),RegionTotalForecastUpperPercentile(11,:),'k-');
        CanterburyLowerPCGraph = plot(FCYear(1:ForecastPeriods),RegionTotalForecastLowerPercentile(11,:),'k-');
        CanterburyMedianGraph = plot(FCYear(1:ForecastPeriods),RegionTotalForecastMedian(11,:),'k-');
        set(CanterburyUpperPCGraph,'LineWidth',3);
        set(CanterburyLowerPCGraph,'LineWidth',3);
        set(CanterburyMedianGraph,'LineWidth',3);
    else
        CanterburyMedian = CanterburyPlotData(1,:);
    end;
    drawnow;
    hold off;
    
    figure(43);
    % Plot of SouthCant
    SouthCantPlotData=zeros(ForecastPeriods,MCRUNS);
    SouthCantPlotData(:,:)=ForecastRegionTotalStore(:,12,:);  
    HistoricalSouthCantGWH = RegionHistAnnualTotals(:,12);
    % and plot
    RegionalResultSpread=plot(FCYear(1:ForecastPeriods),SouthCantPlotData);
    hold on;
    axis([1998 FCYear(ForecastPeriods) 0 2000]);
    title('SouthCant Forecast Growth');
    HistSouthCantPlot = plot(RegionHistYear,HistoricalSouthCantGWH, 'k-');
    set(HistSouthCantPlot,'LineWidth',3);
    if SingleRun~=1;
        SouthCantUpperPCGraph = plot(FCYear(1:ForecastPeriods),RegionTotalForecastUpperPercentile(12,:),'k-');
        SouthCantLowerPCGraph = plot(FCYear(1:ForecastPeriods),RegionTotalForecastLowerPercentile(12,:),'k-');
        SouthCantMedianGraph = plot(FCYear(1:ForecastPeriods),RegionTotalForecastMedian(12,:),'k-');
        set(SouthCantUpperPCGraph,'LineWidth',3);
        set(SouthCantLowerPCGraph,'LineWidth',3);
        set(SouthCantMedianGraph,'LineWidth',3);
    else
        SouthCantMedian = SouthCantPlotData(1,:);
    end;
    drawnow;
    hold off;

    figure(44);
    % Plot of Otago/Southland
    OtagoPlotData=zeros(ForecastPeriods,MCRUNS);
    OtagoPlotData(:,:)=ForecastRegionTotalStore(:,13,:);  
    HistoricalOtagoGWH = RegionHistAnnualTotals(:,13);
    % and plot
    RegionalResultSpread=plot(FCYear(1:ForecastPeriods),OtagoPlotData);
    hold on;
    axis([1998 FCYear(ForecastPeriods) 0 14000]);
    title('Otago/Southland Forecast Growth');
    HistOtagoPlot = plot(RegionHistYear,HistoricalOtagoGWH, 'k-');
    set(HistOtagoPlot,'LineWidth',3);
    if SingleRun~=1;
        OtagoUpperPCGraph = plot(FCYear(1:ForecastPeriods),RegionTotalForecastUpperPercentile(13,:),'k-');
        OtagoLowerPCGraph = plot(FCYear(1:ForecastPeriods),RegionTotalForecastLowerPercentile(13,:),'k-');
        OtagoMedianGraph = plot(FCYear(1:ForecastPeriods),RegionTotalForecastMedian(13,:),'k-');
        set(OtagoUpperPCGraph,'LineWidth',3);
        set(OtagoLowerPCGraph,'LineWidth',3);
        set(OtagoMedianGraph,'LineWidth',3);
    else
        OtagoMedian = OtagoPlotData(1,:);
    end;
    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;
FCRegionGDP = TempFCRegionGDP;

%Warnings back on

warning on MATLAB:divideByZero
warning on MATLAB:nearlySingularMatrix

% Save to GXP forecast .mat file


ForecastYears = FCYear(1:ForecastPeriods);
GXPForecastNames = GXPnames;

if SingleRun == 1
    save 2009GXPForecasts GXPTotalForecastDemand GXPTotalForecastDemandAdj GXPForecastNames ForecastYears Regions NewLoad GXPExistingEmbeddedMW
end;

if (SaveRegionalMCResults == 1 & SingleRun ~= 1)
    save 2009RegionalForecasts Regions FCYear ForecastRegionTotalStore RegionTotalForecastLowerPercentile RegionTotalForecastUpperPercentile RegionTotalForecastMedian
end;

if (SaveGXPMCResults == 1 & SingleRun ~= 1)
    save 2009GXPMCForecasts ForecastGXPTotalStore GXPForecastNames ForecastYears Regions NewLoad GXPExistingEmbeddedMW
end;

toc
disp('Done')