function [paths,model] = getSimulationBegin(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_Begin','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

%% Monetary

% We generate the noise terms and initiate the variables
UMs = unifrnd(0,1,N,T);
momSeries = NaN(N,T+1);
momSeries(:,1) = model.Monetary.regimes(end);

% We generate the scenarios
for dt = 1:T
  momSeries(:,dt+1) = ...
    (momSeries(:,dt) == 1).*( (chains.Monetary.pUU(I) >= UMs(:,dt)).*1   + ...
                              (chains.Monetary.pUU(I) <  UMs(:,dt)).*2 ) + ...
    (momSeries(:,dt) == 2).*( (chains.Monetary.pSU(I) >= UMs(:,dt)).*1   + ...
                              ( (chains.Monetary.pSU(I) < UMs(:,dt)) & (chains.Monetary.pSU(I) + chains.Monetary.pSD(I) >= UMs(:,dt))).*3   + ...
                              (chains.Monetary.pSU(I) + chains.Monetary.pSD(I) < UMs(:,dt)).*2 ) + ...
    (momSeries(:,dt) == 3).*( (chains.Monetary.pDD(I) >= UMs(:,dt)).*3   + ...
                              (chains.Monetary.pDD(I) <  UMs(:,dt)).*2 ) ;
end

paths.Monetary = momSeries;

%% Inflation, Short Rate, Dividends (all together)

% We generate the noise terms and initiate the variables
NIs = normrnd(0,1,N,T);
infSeries = NaN(N,T+1);
infSigma2 = NaN(N,T+1);

chains.Inflation.beta_q = chains.Inflation.alphabetagamma_q - chains.Inflation.alpha_q.*(1+chains.Inflation.gamma_q.^2); chains.Inflation.beta_q = chains.Inflation.beta_q(:);

NRs = normrnd(0,1,N,T);
shrSeries = NaN(N,T+1);
shrSigma2 = NaN(N,T+1);

chains.ShortRate.beta_r = chains.ShortRate.alphabetagamma_r - chains.ShortRate.alpha_r.*(1+chains.ShortRate.gamma_r.^2); chains.ShortRate.beta_r = chains.ShortRate.beta_r(:);

NDs = normrnd(0,1,N,T);
divSeries = NaN(N,T+1);
divSigma2 = NaN(N,T+1);

chains.Dividends.beta_d = chains.Dividends.alphabetagamma_d - chains.Dividends.alpha_d.*(1+chains.Dividends.gamma_d.^2); chains.Dividends.beta_d = chains.Dividends.beta_d(:);

% For each scenario, we fix the parameter value and obtain the past
% conditional variance
for dp = 1:N
  
  % We extract the right parameters and assign them to the Inflation,
  % ShortRate, and Dividend objects
  pvnew = [ ...
      chains.Inflation.mu_q1(I(dp)), ...
      chains.Inflation.mu_q2(I(dp)), ...
      chains.Inflation.mu_q3(I(dp)), ...
      chains.Inflation.a_q(I(dp)), ...
      chains.Inflation.sig2_q(I(dp)), ...
      chains.Inflation.alpha_q(I(dp)), ...
      chains.Inflation.alphabetagamma_q(I(dp)), ...
      chains.Inflation.gamma_q(I(dp)), ...
      chains.Inflation.q0(I(dp)), ...
      chains.Inflation.sig2_q_init(I(dp))];


  [~,z_q,sig2_q] = model.Inflation.getMLEfunction(pvnew);
  eps_q = z_q./sqrt(sig2_q);
  
  pvnew = [ ...
      chains.ShortRate.mu_r1(I(dp)), ...
      chains.ShortRate.mu_r2(I(dp)), ...
      chains.ShortRate.mu_r3(I(dp)), ...
      chains.ShortRate.a_r(I(dp)), ...
      chains.ShortRate.rho_qr(I(dp)), ...
      chains.ShortRate.sig2_r(I(dp)), ...
      chains.ShortRate.alpha_r(I(dp)), ...
      chains.ShortRate.alphabetagamma_r(I(dp)), ...
      chains.ShortRate.gamma_r(I(dp)), ...
      chains.ShortRate.r0(I(dp)), ...
      chains.ShortRate.sig2_r_init(I(dp))];

  [~,z_r,sig2_r] = model.ShortRate.getMLEfunction(pvnew,eps_q);
  eps_r = z_r./sqrt(sig2_r);

  pvnew = [ ...
      chains.Dividends.mu_d1(I(dp)), ...
      chains.Dividends.mu_d2(I(dp)), ...
      chains.Dividends.mu_d3(I(dp)), ...
      chains.Dividends.a_d(I(dp)), ...
      chains.Dividends.rho_qd(I(dp)), ...
      chains.Dividends.rho_rd(I(dp)), ...
      chains.Dividends.sig2_d(I(dp)), ...
      chains.Dividends.alpha_d(I(dp)), ...
      chains.Dividends.alphabetagamma_d(I(dp)), ...
      chains.Dividends.gamma_d(I(dp)), ...
      chains.Dividends.d0(I(dp)), ...
      chains.Dividends.sig2_d_init(I(dp))];

  [~,z_d,sig2_d] = model.Dividends.getMLEfunction(pvnew,eps_q,eps_r,chains.ShortRate.rho_qr(I(dp)));

  % We initialize our series to the last observations
  infSeries(dp,1) = model.Inflation.series(end);
  infSigma2(dp,1) = chains.Inflation.sig2_q(I(dp)) + chains.Inflation.beta_q(I(dp)).*(sig2_q(end) - chains.Inflation.sig2_q(I(dp))) + chains.Inflation.alpha_q(I(dp)).*( (z_q(end) - chains.Inflation.gamma_q(I(dp)).*sqrt(sig2_q(end))).^2 - chains.Inflation.sig2_q(I(dp)).*(1 + chains.Inflation.gamma_q(I(dp)).^2) );

  % We generate the inflation scenarios
  for dt = 1:T
    mu = chains.Inflation.mu_q1(I(dp)).*(momSeries(dp,dt+1) == 1) + ...
         chains.Inflation.mu_q2(I(dp)).*(momSeries(dp,dt+1) == 2) + ...
         chains.Inflation.mu_q3(I(dp)).*(momSeries(dp,dt+1) == 3);
    z  = sqrt(infSigma2(dp,dt)).*NIs(dp,dt);
    infSeries(dp,dt+1) = mu + chains.Inflation.a_q(I(dp)).*(infSeries(dp,dt) - mu) + z;
    infSigma2(dp,dt+1) = boxed(chains.Inflation.sig2_q(I(dp)) + chains.Inflation.beta_q(I(dp)).*(infSigma2(dp,dt) - chains.Inflation.sig2_q(I(dp))) + chains.Inflation.alpha_q(I(dp)).*( (z - chains.Inflation.gamma_q(I(dp)).*sqrt(infSigma2(dp,dt))).^2 - chains.Inflation.sig2_q(I(dp)).*(1+chains.Inflation.gamma_q(I(dp)).^2) ),1e-8,1e-1);
  end
  clear z mu 

  % We generate the inflation scenarios
  shrSeries(dp,1) = model.ShortRate.series(end);
  shrSigma2(dp,1) = chains.ShortRate.sig2_r(I(dp)) + chains.ShortRate.beta_r(I(dp)).*(sig2_r(end) - chains.ShortRate.sig2_r(I(dp))) + chains.ShortRate.alpha_r(I(dp)).*( (z_r(end) - chains.ShortRate.gamma_r(I(dp)).*sqrt(sig2_r(end))).^2 - chains.ShortRate.sig2_r(I(dp)).*(1 + chains.ShortRate.gamma_r(I(dp)).^2) );
  
  % We generate the short rate scenarios
  for dt = 1:T
    mu = chains.ShortRate.mu_r1(I(dp)).*(momSeries(dp,dt+1) == 1) + ...
         chains.ShortRate.mu_r2(I(dp)).*(momSeries(dp,dt+1) == 2) + ...
         chains.ShortRate.mu_r3(I(dp)).*(momSeries(dp,dt+1) == 3);
    z  = sqrt(shrSigma2(dp,dt)).*(chains.ShortRate.rho_qr(I(dp)).*NIs(dp,dt) + sqrt(1-chains.ShortRate.rho_qr(I(dp)).^2).*NRs(dp,dt));
    shrSeries(dp,dt+1) = mu + chains.ShortRate.a_r(I(dp)).*(shrSeries(dp,dt) - mu) + z;
    shrSigma2(dp,dt+1) = boxed(chains.ShortRate.sig2_r(I(dp)) + chains.ShortRate.beta_r(I(dp)).*(shrSigma2(dp,dt) - chains.ShortRate.sig2_r(I(dp))) + chains.ShortRate.alpha_r(I(dp)).*( (z - chains.ShortRate.gamma_r(I(dp)).*sqrt(shrSigma2(dp,dt))).^2 - chains.ShortRate.sig2_r(I(dp)).*(1+chains.ShortRate.gamma_r(I(dp)).^2) ),1e-8,1e-1);
  end
  clear z mu

  % We generate the inflation scenarios
  divSeries(dp,1) = log(model.Dividends.series(end));
  divSigma2(dp,1) = chains.Dividends.sig2_d(I(dp)) + chains.Dividends.beta_d(I(dp)).*(sig2_d(end)-chains.Dividends.sig2_d(I(dp))) + chains.Dividends.alpha_d(I(dp)).*( (z_d(end) - chains.Dividends.gamma_d(I(dp)).*sqrt(sig2_d(end))).^2 - chains.Dividends.sig2_d(I(dp)).*(1 + chains.Dividends.gamma_d(I(dp)).^2) );

  % We generate the dividend scenarios
  C = chol([1,chains.ShortRate.rho_qr(I(dp)),chains.Dividends.rho_qd(I(dp));chains.ShortRate.rho_qr(I(dp)),1,chains.Dividends.rho_rd(I(dp));chains.Dividends.rho_qd(I(dp)),chains.Dividends.rho_rd(I(dp)),1])';
  cs = C(3,:);
  for dt = 1:T
    mu = chains.Dividends.mu_d1(I(dp)).*(momSeries(dp,dt+1) == 1) + ...
         chains.Dividends.mu_d2(I(dp)).*(momSeries(dp,dt+1) == 2) + ...
         chains.Dividends.mu_d3(I(dp)).*(momSeries(dp,dt+1) == 3);
    z  = sqrt(divSigma2(dp,dt)).*(cs(1).*NIs(dp,dt) + cs(2).*NRs(dp,dt) + cs(3).*NDs(dp,dt));
    divSeries(dp,dt+1) = mu + chains.Dividends.a_d(I(dp)).*(divSeries(dp,dt) - mu) + z;
    divSigma2(dp,dt+1) = boxed(chains.Dividends.sig2_d(I(dp)) + chains.Dividends.beta_d(I(dp)).*(divSigma2(dp,dt) - chains.Dividends.sig2_d(I(dp))) + chains.Dividends.alpha_d(I(dp)).*( (z - chains.Dividends.gamma_d(I(dp)).*sqrt(divSigma2(dp,dt))).^2 - chains.Dividends.sig2_d(I(dp)).*(1+chains.Dividends.gamma_d(I(dp)).^2) ),1e-8,1e-1);
  end
  clear cs C mu

end
clear NDs NRs NIs

% 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.Inflation = infSeries;
paths.ShortRate = unshrSeries;
paths.Dividends = exp(divSeries);

%% IndexReturns

% We extract the right parameters and assign them to the IndexReturns 
% objects
pvnew = [ ...
  chains.IndexReturns.mu_y1(I), ...
  chains.IndexReturns.mu_y2(I), ...
  chains.IndexReturns.mu_y3(I), ...
  chains.IndexReturns.sig2_y1(I), ...
  chains.IndexReturns.alpha_y(I), ...
  chains.IndexReturns.alphabetagamma_y(I), ...
  chains.IndexReturns.gamma_y(I), ...
  chains.IndexReturns.sig2_y_init(I)];

% We extract the last innovation and conditional variance
filtered.IndexReturns.z = NaN(N,1);
filtered.IndexReturns.sig2 = NaN(N,1);
for dp = 1:N
  [~,z,sig2] = model.IndexReturns.getMLEfunction(pvnew(dp,:));
  filtered.IndexReturns.z(dp,1) = z(end);
  filtered.IndexReturns.sig2(dp,1) = sig2(end);
end

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

chains.IndexReturns.beta_y = chains.IndexReturns.alphabetagamma_y - chains.IndexReturns.alpha_y.*(1+chains.IndexReturns.gamma_y.^2); chains.IndexReturns.beta_y = chains.IndexReturns.beta_y(:);
retSigma2(:,1) = chains.IndexReturns.sig2_y1(I) + chains.IndexReturns.beta_y(I).*(filtered.IndexReturns.sig2-chains.IndexReturns.sig2_y1(I)) + chains.IndexReturns.alpha_y(I).*( (filtered.IndexReturns.z - chains.IndexReturns.gamma_y(I).*sqrt(filtered.IndexReturns.sig2)).^2 - chains.IndexReturns.sig2_y1(I).*(1 + chains.IndexReturns.gamma_y(I).^2) );

% We generate the scenarios
for dt = 1:T
  mu = chains.IndexReturns.mu_y1(I).*(momSeries(:,dt+1) == 1) + ...
       chains.IndexReturns.mu_y2(I).*(momSeries(:,dt+1) == 2) + ...
       chains.IndexReturns.mu_y3(I).*(momSeries(:,dt+1) == 3);
  z  = sqrt(retSigma2(:,dt)).*NYs(:,dt);
  retSeries(:,dt+1) = mu - retSigma2(:,dt)./2 + z;
  retSigma2(:,dt+1) = boxed(chains.IndexReturns.sig2_y1(I) + chains.IndexReturns.beta_y(I).*(retSigma2(:,dt) - chains.IndexReturns.sig2_y1(I)) + chains.IndexReturns.alpha_y(I).*( (z - chains.IndexReturns.gamma_y(I).*sqrt(retSigma2(:,dt))).^2 - chains.IndexReturns.sig2_y1(I).*(1+chains.IndexReturns.gamma_y(I).^2) ),1e-5,1e0);
end
clear retSigma2 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);
NF2s = normrnd(0,1,N,T);

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

% 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);
  fac2Series(:,dt+1) = chains.ForwardRate.mu_F2(I) + chains.ForwardRate.A_F2(I).*(fac2Series(:,dt) - chains.ForwardRate.mu_F2(I)) + sqrt(chains.ForwardRate.sig2_F2(I)).*NF2s(:,dt);
end
clear NF1s NF2s

% 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)  + chains.ForwardRate.A_f21(I).*fac2Series(:,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)  + chains.ForwardRate.A_f22(I).*fac2Series(:,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)  + chains.ForwardRate.A_f23(I).*fac2Series(:,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)  + chains.ForwardRate.A_f25(I).*fac2Series(:,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)  + chains.ForwardRate.A_f27(I).*fac2Series(:,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) + chains.ForwardRate.A_f210(I).*fac2Series(:,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) + chains.ForwardRate.A_f230(I).*fac2Series(:,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