%% START-------------------------------------------------------------------
close all; clc; clear all;
%% ALLOCATION--------------------------------------------------------------
Eq=5477;
Fi=2533;
Fig=2533*2/3*0.75;
Cash=2533*2/3*0.25;
Fic=2533*(1-2/3);
RE=246;               % unused in the modelling
years=30;             % 30-year period
% years=50;           % 50-year period
nPeriods = years*12;  % number of months
N = 5000;             % number of simulations
sprate = 0.03;        % spending rate
sprate_gdp = 0.07;    % non-oil budget deficit – GDP sprending rate
%% Import data from spreadsheet--------------------------------------------
% Script for importing data from the following spreadsheet:
%
% Workbook: C:\Users\Lenovo\Desktop\Norges\in-the-process\quant\PORTFOLIO.xlsx
%
% Auto-generated by MATLAB on 2019/04/15 23:09:24
[~, ~, raw, dates] = xlsread('C:\Users\Lenovo\Desktop\Norges\in-the-process\quant\PORTFOLIO.xlsx','Eq','A2:AW592','',@convertSpreadsheetExcelDates);
raw(cellfun(@(x) ~isempty(x) && isnumeric(x) && isnan(x),raw)) = {''};
raw = raw(:,[2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49]);
dates = dates(:,1);
R = cellfun(@(x) isempty(x) || (ischar(x) && all(x==' ')),raw);
raw(R) = {0.0}; % Replace blank cells
R = cellfun(@(x) isempty(x) || (ischar(x) && all(x==' ')),dates);
dates(R) = {NaN}; % Replace blank cells
data = reshape([raw{:}],size(raw)); % Create output variable
PORTFOLIOS1 = table; % Create table
% Allocate imported array to column variable names
PORTFOLIOS1.Date = datetime([dates{:,1}].', 'ConvertFrom', 'Excel');
PORTFOLIOS1.Brazil = data(:,1);
PORTFOLIOS1.Chile = data(:,2);
PORTFOLIOS1.Colombia = data(:,3);
PORTFOLIOS1.Canada = data(:,4);
PORTFOLIOS1.Mexico = data(:,5);
PORTFOLIOS1.Peru = data(:,6);
PORTFOLIOS1.USA = data(:,7);
PORTFOLIOS1.Belgium = data(:,8);
PORTFOLIOS1.Finland = data(:,9);
PORTFOLIOS1.France = data(:,10);
PORTFOLIOS1.Greece = data(:,11);
PORTFOLIOS1.Ireland = data(:,12);
PORTFOLIOS1.Italy = data(:,13);
PORTFOLIOS1.Netherlands = data(:,14);
PORTFOLIOS1.Austria = data(:,15);
PORTFOLIOS1.Portugal = data(:,16);
PORTFOLIOS1.Spain = data(:,17);
PORTFOLIOS1.Germany = data(:,18);
PORTFOLIOS1.Euroarea = data(:,19);
PORTFOLIOS1.Denmark = data(:,20);
PORTFOLIOS1.Poland = data(:,21);
PORTFOLIOS1.Russia = data(:,22);
PORTFOLIOS1.UK = data(:,23);
PORTFOLIOS1.Switzerland = data(:,24);
PORTFOLIOS1.Sweden = data(:,25);
PORTFOLIOS1.Czechia = data(:,26);
PORTFOLIOS1.Turkey = data(:,27);
PORTFOLIOS1.Hungary = data(:,28);
PORTFOLIOS1.UAE = data(:,29);
PORTFOLIOS1.Egypt = data(:,30);
PORTFOLIOS1.Israel = data(:,31);
PORTFOLIOS1.Kuwait = data(:,32);
PORTFOLIOS1.Qatar = data(:,33);
PORTFOLIOS1.SouthAfrica = data(:,34);
PORTFOLIOS1.Australia = data(:,35);
PORTFOLIOS1.Philippines = data(:,36);
PORTFOLIOS1.HongKong = data(:,37);
PORTFOLIOS1.India = data(:,38);
PORTFOLIOS1.Indonesia = data(:,39);
PORTFOLIOS1.Japan = data(:,40);
PORTFOLIOS1.China = data(:,41);
PORTFOLIOS1.Malaysia = data(:,42);
PORTFOLIOS1.NewZealand = data(:,43);
PORTFOLIOS1.Pakistan = data(:,44);
PORTFOLIOS1.Singapore = data(:,45);
PORTFOLIOS1.Korea = data(:,46);
PORTFOLIOS1.Taiwan = data(:,47);
PORTFOLIOS1.Thailand = data(:,48);
clearvars data raw dates R; % Clear temporary variables
%--------------------------------------------------------------------------
[~, ~, raw] = xlsread('C:\Users\Lenovo\Desktop\Norges\in-the-process\quant\PORTFOLIO.xlsx','Eqw','A2:AV2');
data = reshape([raw{:}],size(raw)); % Create output variable
PORTFOLIOS2 = table; % Create table
% Allocate imported array to column variable names
PORTFOLIOS2.Brazil = data(:,1);
PORTFOLIOS2.Chile = data(:,2);
PORTFOLIOS2.Colombia = data(:,3);
PORTFOLIOS2.Canada = data(:,4);
PORTFOLIOS2.Mexico = data(:,5);
PORTFOLIOS2.Peru = data(:,6);
PORTFOLIOS2.USA = data(:,7);
PORTFOLIOS2.Belgium = data(:,8);
PORTFOLIOS2.Finland = data(:,9);
PORTFOLIOS2.France = data(:,10);
PORTFOLIOS2.Greece = data(:,11);
PORTFOLIOS2.Ireland = data(:,12);
PORTFOLIOS2.Italy = data(:,13);
PORTFOLIOS2.Netherlands = data(:,14);
PORTFOLIOS2.Austria = data(:,15);
PORTFOLIOS2.Portugal = data(:,16);
PORTFOLIOS2.Spain = data(:,17);
PORTFOLIOS2.Germany = data(:,18);
PORTFOLIOS2.Euroarea = data(:,19);
PORTFOLIOS2.Denmark = data(:,20);
PORTFOLIOS2.Poland = data(:,21);
PORTFOLIOS2.Russia = data(:,22);
PORTFOLIOS2.UK = data(:,23);
PORTFOLIOS2.Switzerland = data(:,24);
PORTFOLIOS2.Sweden = data(:,25);
PORTFOLIOS2.Czechia = data(:,26);
PORTFOLIOS2.Turkey = data(:,27);
PORTFOLIOS2.Hungary = data(:,28);
PORTFOLIOS2.UAE = data(:,29);
PORTFOLIOS2.Egypt = data(:,30);
PORTFOLIOS2.Israel = data(:,31);
PORTFOLIOS2.Kuweit = data(:,32);
PORTFOLIOS2.Qatar = data(:,33);
PORTFOLIOS2.SouthAfrica = data(:,34);
PORTFOLIOS2.Australia = data(:,35);
PORTFOLIOS2.Philippines = data(:,36);
PORTFOLIOS2.HongKong = data(:,37);
PORTFOLIOS2.India = data(:,38);
PORTFOLIOS2.Indonesia = data(:,39);
PORTFOLIOS2.Japan = data(:,40);
PORTFOLIOS2.China = data(:,41);
PORTFOLIOS2.Malaysia = data(:,42);
PORTFOLIOS2.NewZealand = data(:,43);
PORTFOLIOS2.Pakistan = data(:,44);
PORTFOLIOS2.Singapore = data(:,45);
PORTFOLIOS2.Korea = data(:,46);
PORTFOLIOS2.Taiwan = data(:,47);
PORTFOLIOS2.Thailand = data(:,48);
clearvars data raw; % Clear temporary variables
%--------------------------------------------------------------------------
[~, ~, raw, dates] = xlsread('C:\Users\Lenovo\Desktop\Norges\in-the-process\quant\PORTFOLIO.xlsx','Fi','A2:AW592','',@convertSpreadsheetExcelDates);
raw(cellfun(@(x) ~isempty(x) && isnumeric(x) && isnan(x),raw)) = {''};
raw = raw(:,[2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49]);
dates = dates(:,1);
R = cellfun(@(x) isempty(x) || (ischar(x) && all(x==' ')),raw);
raw(R) = {0.0}; % Replace blank cells
R = cellfun(@(x) isempty(x) || (ischar(x) && all(x==' ')),dates);
dates(R) = {NaN}; % Replace blank cells
data = reshape([raw{:}],size(raw)); % Create output variable
PORTFOLIOS3 = table; % Create table
% Allocate imported array to column variable names
PORTFOLIOS3.VarName1 = datetime([dates{:,1}].', 'ConvertFrom', 'Excel');
PORTFOLIOS3.Australia = data(:,1);
PORTFOLIOS3.NewZealand = data(:,2);
PORTFOLIOS3.Brazil = data(:,3);
PORTFOLIOS3.Chile = data(:,4);
PORTFOLIOS3.Mexico = data(:,5);
PORTFOLIOS3.Colombia = data(:,6);
PORTFOLIOS3.USA = data(:,7);
PORTFOLIOS3.Canada = data(:,8);
PORTFOLIOS3.Austria = data(:,9);
PORTFOLIOS3.Belgium = data(:,10);
PORTFOLIOS3.Croatia = data(:,11);
PORTFOLIOS3.Cyprus = data(:,12);
PORTFOLIOS3.Czechia = data(:,13);
PORTFOLIOS3.Denmark = data(:,14);
PORTFOLIOS3.Finland = data(:,15);
PORTFOLIOS3.France = data(:,16);
PORTFOLIOS3.Germany = data(:,17);
PORTFOLIOS3.Hungary = data(:,18);
PORTFOLIOS3.Iceland = data(:,19);
PORTFOLIOS3.Ireland = data(:,20);
PORTFOLIOS3.Italy = data(:,21);
PORTFOLIOS3.Latvia = data(:,22);
PORTFOLIOS3.Lithuania = data(:,23);
PORTFOLIOS3.Luxembourg = data(:,24);
PORTFOLIOS3.Malta = data(:,25);
PORTFOLIOS3.Netherlands = data(:,26);
PORTFOLIOS3.Poland = data(:,27);
PORTFOLIOS3.Portugal = data(:,28);
PORTFOLIOS3.Romania = data(:,29);
PORTFOLIOS3.Russia = data(:,30);
PORTFOLIOS3.Slovakia = data(:,31);
PORTFOLIOS3.Slovenia = data(:,32);
PORTFOLIOS3.Spain = data(:,33);
PORTFOLIOS3.Sweden = data(:,34);
PORTFOLIOS3.Switzerland = data(:,35);
PORTFOLIOS3.Turkey = data(:,36);
PORTFOLIOS3.UK = data(:,37);
PORTFOLIOS3.China = data(:,38);
PORTFOLIOS3.India = data(:,39);
PORTFOLIOS3.Indonesia = data(:,40);
PORTFOLIOS3.Japan = data(:,41);
PORTFOLIOS3.Malaysia = data(:,42);
PORTFOLIOS3.Philippines = data(:,43);
PORTFOLIOS3.Singapore = data(:,44);
PORTFOLIOS3.Korea = data(:,45);
PORTFOLIOS3.Thailand = data(:,46);
PORTFOLIOS3.Israel = data(:,47);
PORTFOLIOS3.SouthAfrica = data(:,48);
clearvars data raw dates R; % Clear temporary variables
%--------------------------------------------------------------------------
[~, ~, raw, dates] = xlsread('C:\Users\Lenovo\Desktop\Norges\in-the-process\quant\PORTFOLIO.xlsx','Fic','A2:C224','',@convertSpreadsheetExcelDates);
raw(cellfun(@(x) ~isempty(x) && isnumeric(x) && isnan(x),raw)) = {''};
raw = raw(:,[3]);
dates = dates(:,1);
R = cellfun(@(x) ~isnumeric(x) && ~islogical(x),raw); % Find non-numeric cells
raw(R) = {NaN}; % Replace non-numeric cells
R = cellfun(@(x) ~isnumeric(x) && ~islogical(x),dates); % Find non-numeric cells
dates(R) = {NaN}; % Replace non-numeric Excel dates with NaN
data = reshape([raw{:}],size(raw)); %Create output variable
PORTFOLIOS4 = table; % Create table
PORTFOLIOS4.date = datetime([dates{:,1}].', 'ConvertFrom', 'Excel');
PORTFOLIOS4.return1 = data(:,1);
clearvars data raw dates R; % Clear temporary variables
%--------------------------------------------------------------------------
[~, ~, raw] = xlsread('C:\Users\Lenovo\Desktop\Norges\in-the-process\quant\PORTFOLIO.xlsx','Fiw','A2:AV2');
data = reshape([raw{:}],size(raw)); %Create output variable
PORTFOLIOS5 = table; % Create table
% Allocate imported array to column variable names
PORTFOLIOS5.Australia = data(:,1);
PORTFOLIOS5.NewZealand = data(:,2);
PORTFOLIOS5.Brazil = data(:,3);
PORTFOLIOS5.Chile = data(:,4);
PORTFOLIOS5.Mexico = data(:,5);
PORTFOLIOS5.Colombia = data(:,6);
PORTFOLIOS5.USA = data(:,7);
PORTFOLIOS5.Canada = data(:,8);
PORTFOLIOS5.Austria = data(:,9);
PORTFOLIOS5.Belgium = data(:,10);
PORTFOLIOS5.Croatia = data(:,11);
PORTFOLIOS5.Cyprus = data(:,12);
PORTFOLIOS5.Czechia = data(:,13);
PORTFOLIOS5.Denmark = data(:,14);
PORTFOLIOS5.Finland = data(:,15);
PORTFOLIOS5.France = data(:,16);
PORTFOLIOS5.Germany = data(:,17);
PORTFOLIOS5.Hungary = data(:,18);
PORTFOLIOS5.Iceland = data(:,19);
PORTFOLIOS5.Ireland = data(:,20);
PORTFOLIOS5.Italy = data(:,21);
PORTFOLIOS5.Latvia = data(:,22);
PORTFOLIOS5.Lithuania = data(:,23);
PORTFOLIOS5.Luxembourg = data(:,24);
PORTFOLIOS5.Malta = data(:,25);
PORTFOLIOS5.Netherlands = data(:,26);
PORTFOLIOS5.Poland = data(:,27);
PORTFOLIOS5.Portugal = data(:,28);
PORTFOLIOS5.Romania = data(:,29);
PORTFOLIOS5.Russia = data(:,30);
PORTFOLIOS5.Slovakia = data(:,31);
PORTFOLIOS5.Slovenia = data(:,32);
PORTFOLIOS5.Spain = data(:,33);
PORTFOLIOS5.Sweden = data(:,34);
PORTFOLIOS5.Switzerland = data(:,35);
PORTFOLIOS5.Turkey = data(:,36);
PORTFOLIOS5.UK = data(:,37);
PORTFOLIOS5.China = data(:,38);
PORTFOLIOS5.India = data(:,39);
PORTFOLIOS5.Indonesia = data(:,40);
PORTFOLIOS5.Japan = data(:,41);
PORTFOLIOS5.Malaysia = data(:,42);
PORTFOLIOS5.Philippines = data(:,43);
PORTFOLIOS5.Singapore = data(:,44);
PORTFOLIOS5.Korea = data(:,45);
PORTFOLIOS5.Thailand = data(:,46);
PORTFOLIOS5.Israel = data(:,47);
PORTFOLIOS5.SouthAfrica = data(:,48);
clearvars data raw; % Clear temporary variables
%--------------------------------------------------------------------------
[~, ~, raw, dates] = xlsread('C:\Users\Lenovo\Desktop\Norges\in-the-process\quant\PORTFOLIO.xlsx','gdp_infl','A2:D50','',@convertSpreadsheetExcelDates);
raw(cellfun(@(x) ~isempty(x) && isnumeric(x) && isnan(x),raw)) = {''};
raw = raw(:,[2,3,4]);
dates = dates(:,1);
R = cellfun(@(x) ~isnumeric(x) && ~islogical(x),raw); % Find non-numeric cells
raw(R) = {NaN}; % Replace non-numeric cells
R = cellfun(@(x) ~isnumeric(x) && ~islogical(x),dates); % Find non-numeric cells
dates(R) = {NaN}; % Replace non-numeric Excel dates with NaN
data = reshape([raw{:}],size(raw)); % Create output variable
PORTFOLIOS6 = table; % Create table
% Allocate imported array to column variable names
PORTFOLIOS6.VarName1 = datetime([dates{:,1}].', 'ConvertFrom', 'Excel');
PORTFOLIOS6.gdp_growth = data(:,1);
PORTFOLIOS6.ann_defl = data(:,2);
PORTFOLIOS6.ann_cpi = data(:,3);
clearvars data raw dates R; % Clear temporary variables
%% Arrayisation------------------------------------------------------------

Eq_returns=table2array(PORTFOLIOS1(:,2:end));
Eq_weights=table2array(PORTFOLIOS2);
Fig_returns=table2array(PORTFOLIOS3(:,2:end));
Fic_returns=1+table2array(PORTFOLIOS4(:,2:end));
Fig_weights=table2array(PORTFOLIOS5);
Gdp_infl=table2array(PORTFOLIOS6(:,2:end));

Portfolio=horzcat(1+Eq_returns, 1+Fig_returns);
%% Set up------------------------------------------------------------------
Eq_abs=Eq;
Fig_abs=Fig+Cash;
Fic_abs=Fic;
Abs_Fic=Fic;
Eq_abs_post=Eq;
Fig_abs_post=Fig;
Fic_abs_post=Fic; 
Annual_portfolio=Eq_abs+Fig_abs+Fic_abs;
Eq_ovrl_weights=(1/sum(Eq_weights)*0.7).*Eq_weights(1,:);   % target Equity allocation
Fi_ovrl_weights=(0.75/sum(Fig_weights)*0.2).*Fig_weights(1,:); % target Fixed Income allocation
weights=[Eq_ovrl_weights,Fi_ovrl_weights];
Abs_weights=Annual_portfolio.*weights;
Ongoing_Annual_portfolio=(Eq+Fi)*(ones(1,N));
withdrawal=0.03*(Ongoing_Annual_portfolio).*(ones(1,N));
portfolio_returns=1.*(ones(1,N));
GDP=3450;
Infl=1;
kk=randperm(590)';
%% Simulation--------------------------------------------------------------
for u = 1:N
    Annual_portfolio=Eq+Fig+Cash+Fic;
    Abs_weights=Annual_portfolio.*weights;
    Abs_Fic=2533*1/3;
    Cash=2533*2/3*(1-0.75);
    for i=1:nPeriods
        k=ceil(rand*591);
        l=ceil(rand*223);
        Portfolio(591+i,:)= Portfolio(k,:);     % with replacement
%         Portfolio(591+i,:)= Portfolio(kk(i),:); % without replacement
        Fic_returns(223+i,1)= Fic_returns(l,1);
    end
    for i=1:years
        Gdp_infl(49+i,1)=normrnd(1.5,0.75); %GDP growth 
        rho = 0.3;
        Gdp_infl(49+i,2) = rho*Gdp_infl(49+i,1)+((1-rho^2)^0.5)*normrnd(0,1); %inflation rate while preserving correlation
        if Gdp_infl(49+i,2) >= 2.5
            Gdp_infl(49+i,3) = 2.5;
        else
            Gdp_infl(49+i,3)=Gdp_infl(49+i,2);
        end
    end
    for i=1:years
        for j=1:12
            Abs_weights(1+j,:)=(Portfolio(591+(i-1)*12+j,1:96)).*Abs_weights(j,:);
            Abs_Fic(1+j,1)=(Fic_returns(223+(i-1)*12+j))*Abs_Fic(j,1);
        end

        GDP(1+i,1)=GDP(i,1)*(1+Gdp_infl(49+i,1)/100); % GDP forecast
        Infl(1+i,1)=Infl(i,1)*(1+(Gdp_infl(49+i,3))/100); % Inflation forecast
        Cash(1+i,1)=Cash(i,1)*(Infl(1+i,1)/Infl(i,1));        
        
        Eq_abs(1+i,1)=sum(Abs_weights(1+j,1:48));% Equity
        Fig_abs(1+i,1)=sum(Abs_weights(1+j,49:96))+Cash(1+i,1); % FI Government
        Fic_abs(1+i,1)=Abs_Fic(1+j,1); % FI Corporate
        
        Annual_portfolio(1+i,1)=Eq_abs(1+i,1)+Fig_abs(1+i,1)+Fic_abs(1+i,1); % All
        %------------Constant----------------------------------------------
%           withdrawal(1+i,u)=Annual_portfolio(1+i,1)*sprate; %(3% Rule)
        %------------Average-----------------------------------------------
        %-----------------------3-year-------------------------------------
%         if i == 1
%             withdrawal(1+i,u)=(Annual_portfolio(1+i,1)+Annual_portfolio(i,1))*sprate/2;
%         else
%             withdrawal(1+i,u)=(Annual_portfolio(1+i,1)+Annual_portfolio(i,1)+...
%             Annual_portfolio(i-1,1))*sprate/3;
%         end
%-------------------------------5-year-------------------------------------        
%         if i == 1
%             withdrawal(1+i,u)=(Annual_portfolio(1+i,1)+Annual_portfolio(i,1))*sprate/2;
%         elseif i == 2
%             withdrawal(1+i,u)=(Annual_portfolio(1+i,1)+Annual_portfolio(i,1)+...
%             Annual_portfolio(i-1,1))*sprate/3;
%         elseif i == 3
%             withdrawal(1+i,u)=(Annual_portfolio(1+i,1)+Annual_portfolio(i,1)+...
%             Annual_portfolio(i-1,1)+Annual_portfolio(i-2,1))*sprate/4;
%         else
%             withdrawal(1+i,u)=(Annual_portfolio(1+i,1)+Annual_portfolio(i,1)+...
%             Annual_portfolio(i-1,1)+Annual_portfolio(i-2,1)+Annual_portfolio(i-3,1))*sprate/5;
%         end
        %--------------Ratchet---------------------------------------------
%         if i==1
%             withdrawal(1+i,u)=max(237,Annual_portfolio(1+i,1)*sprate);
%         else
%             withdrawal(1+i,u)=max(withdrawal(i,u),Annual_portfolio(1+i,1)*sprate);
%         end
%         --------------GDP rule--------------------------------------------
        withdrawal(1+i,u)=GDP(1+i,1)*sprate_gdp;
        %------------------------------------------------------------------
        Annual_portfolio(1+i,1)=Annual_portfolio(1+i,1)-withdrawal(1+i,u); % after withdrawal
        old_abs_weights(1,:)=Abs_weights(1+j,:).*(1-sprate); 
        new_abs_weights=Annual_portfolio(1+i,1).*weights(1,:);
        reb_costs=0.001*sum(abs(old_abs_weights-new_abs_weights)); % rebalancing costs
        man_costs=0.0006*Annual_portfolio(1+i,1); % management costs
        Annual_portfolio(1+i,1)=Annual_portfolio(1+i,1)-reb_costs-man_costs; %after rebalancing costs and management costs
        Ongoing_Annual_portfolio(1+i,u)=Annual_portfolio(1+i,1); % closing portfolio value
        Abs_weights=Annual_portfolio(1+i,1)*weights(1,:); %new allocation
        Abs_Fic=Annual_portfolio(1+i,1)*0.1; % new allocation to corporate bonds
        Cash(1+i,1)=Annual_portfolio(1+i,1)*0.05;
        Eq_abs_post(1+i,1)=sum(Abs_weights(1,1:48)); % rebalanced Equity weight
        Fig_abs_post(1+i,1)=sum(Abs_weights(1,49:96)); % rebalanced FI Government weight
        Fic_abs_post(1+i,1)=Abs_Fic; % rebalanced FI Corporate weight
    end 
end

for k=1:years+1
    Eq_share(k,1)=Eq_abs(k,1)/Annual_portfolio(k,1);
    Fig_share(k,1)=Fig_abs(k,1)/Annual_portfolio(k,1);
    Fic_share(k,1)=Fic_abs(k,1)/Annual_portfolio(k,1);
    Eq_share_post(k,1)=Eq_abs_post(k,1)/Annual_portfolio(k,1);
    Fig_share_post(k,1)=(Fig_abs_post(k,1)+Cash(k,1))/Annual_portfolio(k,1);
    Fic_share_post(k,1)=Fic_abs_post(k,1)/Annual_portfolio(k,1);
end
for u =1:N
    for m=1:years
        portfolio_returns(m,u)=Ongoing_Annual_portfolio(1+m,u)/Ongoing_Annual_portfolio(m,u)-1;
    end
end
ll=mean(mean(portfolio_returns)');
real_Ongoing_Annual_portfolio(:,:)=Ongoing_Annual_portfolio(:,:)./Infl(:,1);
real_withdrawal(:,:)=withdrawal(:,:)./Infl(:,1);
share=[Eq_share,Fig_share,Fic_share,Eq_share_post,Fig_share_post,Fic_share_post];
% Portfolio value
naverage=mean(Ongoing_Annual_portfolio,2);
nperc_10=prctile(Ongoing_Annual_portfolio,10,2);
nperc_25=prctile(Ongoing_Annual_portfolio,25,2);
nperc_50=prctile(Ongoing_Annual_portfolio,50,2);
nperc_75=prctile(Ongoing_Annual_portfolio,75,2);
nperc_90=prctile(Ongoing_Annual_portfolio,90,2);
ndescr=[nperc_10,nperc_25,nperc_50,naverage,nperc_75,nperc_90];

raverage=mean(real_Ongoing_Annual_portfolio,2);
rperc_10=prctile(real_Ongoing_Annual_portfolio,10,2);
rperc_25=prctile(real_Ongoing_Annual_portfolio,25,2);
rperc_50=prctile(real_Ongoing_Annual_portfolio,50,2);
rperc_75=prctile(real_Ongoing_Annual_portfolio,75,2);
rperc_90=prctile(real_Ongoing_Annual_portfolio,90,2);
rdescr=[rperc_10,rperc_25,rperc_50,raverage,rperc_75,rperc_90];

paverage=mean(real_withdrawal,2);
pperc_10=prctile(real_withdrawal,10,2);
pperc_25=prctile(real_withdrawal,25,2);
pperc_50=prctile(real_withdrawal,50,2);
pperc_75=prctile(real_withdrawal,75,2);
pperc_90=prctile(real_withdrawal,90,2);
pdescr=[pperc_10,pperc_25,pperc_50,paverage,pperc_75,pperc_90];
%% Sustainability----------------------------------------------------------
%----Criterion 1-------------------------------------------------------
fund_max=Ongoing_Annual_portfolio(1,:);
ss1_i25=0*ones(1,N);
ss1_i50=0*ones(1,N);
ss1_m25=0*ones(1,N);
ss1_m50=0*ones(1,N);

for u = 1:N
    for i=1:years
        if real_Ongoing_Annual_portfolio(1+i,u) > fund_max(1,u)
            fund_max(1,u)=real_Ongoing_Annual_portfolio(1+i,u);
        end
        %------------------------------------------------------------------
        if real_Ongoing_Annual_portfolio(years+1,u) < 0.75*Ongoing_Annual_portfolio(1,u)
            ss1_i25(1+i,u)=1;
        else
            ss1_i25(1+i,u)=0;
        end
        if real_Ongoing_Annual_portfolio(years+1,u) < 0.5*Ongoing_Annual_portfolio(1,u)
            ss1_i50(1+i,u)=1;
        else
            ss1_i50(1+i,u)=0;
        end
        %------------------------------------------------------------------
        if real_Ongoing_Annual_portfolio(years+1,u) < 0.75*fund_max(1,u)
            ss1_m25(1+i,u)=1;
        else
            ss1_m25(1+i,u)=0;
        end
        if real_Ongoing_Annual_portfolio(years+1,u) < 0.5*fund_max(1,u)
            ss1_m50(1+i,u)=1;
        else
            ss1_m50(1+i,u)=0;
        end
    end
end
dur_i25=mean(sum(ss1_i25==0));
dur_i50=mean(sum(ss1_i50==0));

crit1_i25=1-sum(ss1_i25(years+1,:)==0)/N;
crit1_i50=1-sum(ss1_i50(years+1,:)==0)/N;
crit1_m25=1-sum(ss1_m25(years+1,:)==0)/N;
crit1_m50=1-sum(ss1_m50(years+1,:)==0)/N;
% check
crit1_i25=sum(ss1_i25(years+1,:)==1)/N;
crit1_i50=sum(ss1_i50(years+1,:)==1)/N;
crit1_m25=sum(ss1_m25(years+1,:)==1)/N;
crit1_m50=sum(ss1_m50(years+1,:)==1)/N;
%----Criterion 2-------------------------------------------------------
ss2=0*ones(1,N);
for u = 1:N
        if mean(real_withdrawal(:,u))<237
            ss2(2,u)=1;
        else
            ss2(2,u)=0;
        end
end
crit2=sum(ss2(2,:)==1)/N;
%----Criterion 3-----------------------------------------------------------
crit3_p10=std(real_withdrawal(11,:))/mean(real_withdrawal(11,:));
crit3_p20=std(real_withdrawal(21,:))/mean(real_withdrawal(21,:));
crit3_p30=std(real_withdrawal(31,:))/mean(real_withdrawal(31,:));

for i=1:years
    crit3_p(i,1)=std(real_withdrawal(i+1,:))/mean(real_withdrawal(i+1,:));
end
real_withdrawal_change=0*ones(1,N);
ss3=0*ones(1,N);
for u = 1:N
    for i=1:years
        real_withdrawal_change(i,u) = log(real_withdrawal(1+i,u)/real_withdrawal(i,u));
    end
    ss3(2,u)=std(real_withdrawal_change(:,u));
end

crit3_ch=mean(ss3(2,:)); 

%distribution of crit3

crit3_perc_10=prctile(ss3,10,2);
crit3_perc_25=prctile(ss3,25,2);
crit3_perc_50=prctile(ss3,50,2);
crit3_mean=mean(ss3,2);
crit3_perc_75=prctile(ss3,75,2);
crit3_perc_90=prctile(ss3,90,2);
crit3_descr=[crit3_perc_10,crit3_perc_25,crit3_perc_50,crit3_mean,crit3_perc_75,crit3_perc_90];
%% Visual------------------------------------------------------------------
m1 = median(real_Ongoing_Annual_portfolio(31,:));
m2 = median(real_withdrawal(31,:));
% Correlation between GDP growth and Inflation
[n1,ctr1] = hist(Gdp_infl(:,1),10);
[n2,ctr2] = hist(Gdp_infl(:,3),10);
subplot(2,2,2);
plot(Gdp_infl(:,3), Gdp_infl(:,1),'.');
axis([-0.5 3.0 -0.5 3.0]);
h1 = gca;
xlabel('Inflation CPI');
ylabel('GDP growth');
subplot(2,2,4);
bar(ctr1,-n1,1);
axis([-0.5 3.0 -max(n1)*1.1 0]);
axis('off');
h2 = gca;
subplot(2,2,1);
barh(ctr2,-n2,1);
axis([-max(n2)*1.1 0 -0.5 3.0]);
axis('off');
h3 = gca;
h1.Position = [0.35 0.35 0.55 0.55];
h2.Position = [.35 .1 .55 .15];
h3.Position = [.1 .35 .15 .55];
colormap([.8 .8 1]);
ss=corr(Gdp_infl(:,1),Gdp_infl(:,3))
