function [paths,model] = getSimulationExtWilkie(T,N)
% This function generates scenarios based on the estimated Markov chains
% INPUTS:   T               Horizon of simulation (in months)
%           N               Number of paths
%
% OUTPUTS:  paths           Paths
%           model           Model object

load('results/chains_ExtWilkie','chains','model'); 

% If we want more than 100,000 paths, we need to sample with replacement
if N < 100000
  I = randsample(1:100000,N,false);
else
  I = randsample(1:100000,N,true);
end

%% Inflation

% We generate the noise terms and initiate the variables
NIs = normrnd(0,1,N,T);
infSeries = NaN(N,T+1);
infSeries(:,1) = model.Inflation.series(end);

% We generate the scenarios
for dt = 1:T
  mu = chains.Inflation.mu_q1(I);
  z  = sqrt(chains.Inflation.sig2_q(I)).*NIs(:,dt);
  infSeries(:,dt+1) = mu + chains.Inflation.a_q(I).*(infSeries(:,dt) - mu) + z;
end
clear z mu 

paths.Inflation = infSeries;
clear infSeries

%% ShortRate

% We generate the noise terms and initiate the variables
NRs = normrnd(0,1,N,T);
shrSeries = NaN(N,T+1);
shrSeries(:,1) = model.ShortRate.series(end);

% We generate the scenarios
for dt = 1:T
  mu = chains.ShortRate.mu_r1(I);
  z  = sqrt(chains.ShortRate.sig2_r(I)).*(chains.ShortRate.rho_qr(I).*NIs(:,dt) + sqrt(1-chains.ShortRate.rho_qr(I).^2).*NRs(:,dt));
  shrSeries(:,dt+1) = mu + chains.ShortRate.a_r(I).*(shrSeries(:,dt) - mu) + z;
end
clear z mu

% We transformed our rates (based on the inverse transformation)
a0 = model.ShortRate.rbar - (model.ShortRate.rbar - model.ShortRate.c)*log(model.ShortRate.rbar - model.ShortRate.c);
a1 = model.ShortRate.rbar - model.ShortRate.c;
unshrSeries = shrSeries.*(shrSeries >= model.ShortRate.rbar) + (exp((shrSeries - a0)./a1) + model.ShortRate.c).*(shrSeries < model.ShortRate.rbar);

paths.ShortRate = unshrSeries;

%% Dividends

% We generate the noise terms and initiate the variables
NDs = normrnd(0,1,N,T);
divSeries = NaN(N,T+1);
divSeries(:,1) = log(model.Dividends.series(end));

% We generate the scenarios
for df = 1:N
  C = chol([1,chains.ShortRate.rho_qr(I(df)),chains.Dividends.rho_qd(I(df));chains.ShortRate.rho_qr(I(df)),1,chains.Dividends.rho_rd(I(df));chains.Dividends.rho_qd(I(df)),chains.Dividends.rho_rd(I(df)),1])';
  cs(df,:) = C(3,:);
end

for dt = 1:T
  mu = chains.Dividends.mu_d1(I);
  z  = sqrt(chains.Dividends.sig2_d(I)).*(cs(:,1).*NIs(:,dt) + cs(:,2).*NRs(:,dt) + cs(:,3).*NDs(:,dt));
  divSeries(:,dt+1) = mu + chains.Dividends.a_d(I).*(divSeries(:,dt) - mu) + z;
end
clear cs C mu divSigma2 NDs NRs NIs

paths.Dividends = exp(divSeries);
clear divSeries

%% IndexReturns

% We generate the noise terms and initiate the variables
NYs = normrnd(0,1,N,T);
retSeries = NaN(N,T+1);
retSeries(:,1) = model.IndexReturns.series(end);

% We generate the scenarios
for dt = 1:T
  mu = chains.IndexReturns.mu_y1(I);
  z  = sqrt(chains.IndexReturns.sig2_y1(I)).*NYs(:,dt);
  retSeries(:,dt+1) = mu - chains.IndexReturns.sig2_y1(I)./2 + z;
end
clear NYs mu z

paths.IndexReturns = unshrSeries./12 + retSeries;
clear retSeries

%% ForwardRate

% We generate the noise terms and initiate the variables
NF1s = normrnd(0,1,N,T);

% We obtain the factors
y  = model.ForwardRate.transfseries;
fac1Series = NaN(N,T+1);
fac1Series(:,1) = y(end,7) - y(end,1);

% We generate the factors
for dt = 1:T
  fac1Series(:,dt+1) = chains.ForwardRate.mu_F1(I) + chains.ForwardRate.A_F1(I).*(fac1Series(:,dt) - chains.ForwardRate.mu_F1(I)) + sqrt(chains.ForwardRate.sig2_F1(I)).*NF1s(:,dt);
end
clear NF1s 

% We generate the noise terms and initiate the variables
Nf1s  = normrnd(0,1,N,T);
Nf2s  = normrnd(0,1,N,T);
Nf3s  = normrnd(0,1,N,T);
Nf5s  = normrnd(0,1,N,T);
Nf7s  = normrnd(0,1,N,T);
Nf10s = normrnd(0,1,N,T);
Nf30s = normrnd(0,1,N,T);

f1Series  = NaN(N,T+1);
f2Series  = NaN(N,T+1);
f3Series  = NaN(N,T+1);
f5Series  = NaN(N,T+1);
f7Series  = NaN(N,T+1);
f10Series = NaN(N,T+1);
f30Series = NaN(N,T+1);

f1Series(:,1)  = y(end,1);
f2Series(:,1)  = y(end,2);
f3Series(:,1)  = y(end,3);
f5Series(:,1)  = y(end,4);
f7Series(:,1)  = y(end,5);
f10Series(:,1) = y(end,6);
f30Series(:,1) = y(end,7);

% We generate the rates
for dt = 1:T
  f1Series(:,dt+1)  = chains.ForwardRate.mu_f1(I)  + chains.ForwardRate.A_f11(I).*fac1Series(:,dt+1)   + sqrt(chains.ForwardRate.sig2_f1(I)).*Nf1s(:,dt);
  f2Series(:,dt+1)  = chains.ForwardRate.mu_f2(I)  + chains.ForwardRate.A_f12(I).*fac1Series(:,dt+1)   + sqrt(chains.ForwardRate.sig2_f2(I)).*Nf2s(:,dt);
  f3Series(:,dt+1)  = chains.ForwardRate.mu_f3(I)  + chains.ForwardRate.A_f13(I).*fac1Series(:,dt+1)   + sqrt(chains.ForwardRate.sig2_f3(I)).*Nf3s(:,dt);
  f5Series(:,dt+1)  = chains.ForwardRate.mu_f5(I)  + chains.ForwardRate.A_f15(I).*fac1Series(:,dt+1)   + sqrt(chains.ForwardRate.sig2_f5(I)).*Nf5s(:,dt);
  f7Series(:,dt+1)  = chains.ForwardRate.mu_f7(I)  + chains.ForwardRate.A_f17(I).*fac1Series(:,dt+1)   + sqrt(chains.ForwardRate.sig2_f7(I)).*Nf7s(:,dt);
  f10Series(:,dt+1) = chains.ForwardRate.mu_f10(I) + chains.ForwardRate.A_f110(I).*fac1Series(:,dt+1)  + sqrt(chains.ForwardRate.sig2_f10(I)).*Nf10s(:,dt);
  f30Series(:,dt+1) = chains.ForwardRate.mu_f30(I) + chains.ForwardRate.A_f130(I).*fac1Series(:,dt+1)  + sqrt(chains.ForwardRate.sig2_f30(I)).*Nf30s(:,dt);
end
clear Nf1s Nf2s Nf3s Nf5s Nf7s Nf10s Nf30s

% We transformed our rates (based on the inverse transformation)
a0 = model.ForwardRate.fbar - (model.ForwardRate.fbar - model.ForwardRate.c)*log(model.ForwardRate.fbar - model.ForwardRate.c);
a1 = model.ForwardRate.fbar - model.ForwardRate.c;
unf1Series  = f1Series + shrSeries;
unf1Series  = unf1Series.*(unf1Series >= model.ForwardRate.fbar) + (exp((unf1Series - a0)./a1) + model.ForwardRate.c).*(unf1Series < model.ForwardRate.fbar);
unf2Series  = f2Series + shrSeries;
unf2Series  = unf2Series.*(unf2Series >= model.ForwardRate.fbar) + (exp((unf2Series - a0)./a1) + model.ForwardRate.c).*(unf2Series < model.ForwardRate.fbar);
unf3Series  = f3Series + shrSeries;
unf3Series  = unf3Series.*(unf3Series >= model.ForwardRate.fbar) + (exp((unf3Series - a0)./a1) + model.ForwardRate.c).*(unf3Series < model.ForwardRate.fbar);
unf5Series  = f5Series + shrSeries;
unf5Series  = unf5Series.*(unf5Series >= model.ForwardRate.fbar) + (exp((unf5Series - a0)./a1) + model.ForwardRate.c).*(unf5Series < model.ForwardRate.fbar);
unf7Series  = f7Series + shrSeries;
unf7Series  = unf7Series.*(unf7Series >= model.ForwardRate.fbar) + (exp((unf7Series - a0)./a1) + model.ForwardRate.c).*(unf7Series < model.ForwardRate.fbar);
unf10Series = f10Series + shrSeries;
unf10Series = unf10Series.*(unf10Series >= model.ForwardRate.fbar) + (exp((unf10Series - a0)./a1) + model.ForwardRate.c).*(unf10Series < model.ForwardRate.fbar);
unf30Series = f30Series + shrSeries;
unf30Series = unf30Series.*(unf30Series >= model.ForwardRate.fbar) + (exp((unf30Series - a0)./a1) + model.ForwardRate.c).*(unf30Series < model.ForwardRate.fbar);

r0  =  unshrSeries;
r1  = (unf1Series.*0.75 + r0/4).*1;
r2  = (unf2Series + r1.*1)./2;
r3  = (unf3Series + r2.*2)./3;
r5  = (unf5Series.*2 + r3.*3)./5;
r7  = (unf7Series.*2 + r5.*5)./7;
r10 = (unf10Series.*3 + r7.*7)./10;
r30 = (unf30Series.*20 + r10.*10)./30;
clear unf1Series unf2Series unf3Series unf5Series unf7Series unf10Series unf30Series

paths.ForwardRate.r1  = r1;
paths.ForwardRate.r2  = r2;
paths.ForwardRate.r3  = r3;
paths.ForwardRate.r5  = r5;
paths.ForwardRate.r7  = r7;
paths.ForwardRate.r10 = r10;
paths.ForwardRate.r30 = r30;

end