function y = y_sim(X,PhiQ,kinfQ,Omega,sige2,rlb,mats,dt)
% simulate forward rates and yields from our shadow rate DTSM 
% with drifting trends
% At each time t, given state variables and parameters, we compute 
% estimated forward rate using eq(7) in WX (2016)  
% f_{n,t} = f_{n,t,estimated} + innovation
% According to eq(6) in WX, f{n,t} = (n+1)*y_{n+1,t} - n*y_{n,t}
% Hence, f_{0,t} = y_{1,t} = s_t = \delta0 + \delta1 * X_t
% y_{2,t} = 1/2 * [y_{1,t} + f{1,t}]
% y_{3,t} = 1/3 * [2*y{2,t}+f{2,t}] = 1/3* [y_{1,t}+f{1,t}+f_{2,t}]
%...
% y_{n,t} = 1/n *[y_{1,t} + f_{1,t} + f_{2,t} +...+ f{n-1,t}]

% INPUT:
%       X:      3xT spanned state variables determining forward rates
%       \Phi^Q: 3x3 diaganol matrix 
%       kinfQ:  scalar, determine \mu^Q
%       \Omega: 3x3 matrix
%       sige2:  scalar, variance of innovation
%       rlb:    scalar, zero lower bound
%       mats:   Kx1, 
%       dt:     scalar
%
% OUTPUT:
%        yields: jjxT, simulated yields
%        fwdr: (jj+1)xT, simulated forward rates


T = size(X,2);
% K = length(mats);
J = mats/dt ;
JJ = 1:1:max(J) ;
loads = SRTSM_loadings(kinfQ, PhiQ, Omega, JJ', dt);
aJ = loads.an;            % jjx1
bJ = loads.bn;            % jjx3
sigmasJ = loads.SigmaQ;   % jjx1
fwdr = nan(length(JJ)+1,T);
innov = sqrt(sige2)*randn(length(JJ),T);

sr = 0 + ones(1,3)*X ;  % 1xT short rate at each point
fwdr(1,:) = sr;
musJ = aJ + bJ*X;  % jjxT
sigmasJ_mat = repmat(sigmasJ,1,T);  % jjxT
z1_temp = (musJ-rlb)./sigmasJ_mat;
z2_temp = rlb + (musJ-rlb).*normcdf(z1_temp) + sigmasJ_mat.*normpdf(z1_temp);
fwdr(2:end,:) = z2_temp + innov;  % annulaized forward rate
% ensure fwdr no less than zero
fwdr(fwdr<0)=0;
y_temp = cumsum(fwdr,1);
dt_mat = repmat(JJ',1,T);
y_mat = y_temp(1:end-1,:)./dt_mat;  % annualized yields

y.yields = y_mat;
y.fwdr = fwdr;


end
