function y = SRTSM_loadings(kinfQ, PhiQ, Sigma, mats, dt)
% generate loadings for shadow rate in Wu and Xia (2016 JMCB)
% mats are measure in month
% dt is the length of each period, for quarterly data, it is 3
% output is the loadings of the shadow yields 
% yields are express in a per period basis

% delta1 = ones(3,1) ;
% muQ = zeros(3,1) ;
% muQ(1,1) = kinfQ ;
J = mats/dt ;
% m = length(J) ;
% N = size(K1QX,1) ;
% Phi = K1QX + eye(N) ;
lambd1 = PhiQ(1,1) ;
lambd2 = PhiQ(2,2) ;
lambd3 = PhiQ(3,3) ;
sigma11 = Sigma(1,1); 
sigma22 = Sigma(2,2); 
sigma33 = Sigma(3,3); 
sigma12 = Sigma(1,2); 
sigma13 = Sigma(1,3); 
sigma23 = Sigma(2,3);
JJ = 1:1:max(J) ;   % max(J) = K
JJ = JJ' ;    % Kx1
% Compute aJ, bJ and sigmasJ
bJ = [lambd1.^J lambd2.^J lambd3.^J] ;  % Jx3
anbar = kinfQ * lambd1.^(JJ-1) ;  
anbar = cumsum(anbar,1) ;  % Kx1
Bn = [lambd1.^(JJ-1) lambd2.^(JJ-1) lambd3.^(JJ-1)]; % Kx3
Bn = cumsum(Bn,1);
zt = [lambd1.^(2*(JJ-1)) lambd2.^(2*(JJ-1)) lambd3.^(2*(JJ-1)) 2*(lambd1*lambd2).^(JJ-1) 2*(lambd1*lambd3).^(JJ-1) 2*(lambd2*lambd3).^(JJ-1) ] ;  % Kx6
zt = cumsum(zt,1) ;
sigmas2 = zt * [sigma11; sigma22; sigma33; sigma12; sigma13; sigma23] ;
aJ = NaN(length(J),1);
sigmaJ2 = NaN(length(J),1);
for i = 1:1:length(J)
    % aJ(i) = anbar(J(i)) - 0.5 * Bn(J(i),:) * Sigma * Bn(J(i),:)';
    aJ(i) = anbar(J(i)) - 0.5 * Bn(J(i),:) * Sigma * Bn(J(i),:)'/1200 ;
%     aJ(i) = anbar(J(i)) - 0.5 * Bn(J(i),:) * Sigma * Bn(J(i),:)'/400 ;
    % aJ(i) = delta0 - 0.5 * Bn(J(i),:) * Sigma * Bn(J(i),:)'/400 ;
    sigmaJ2(i) = sigmas2(J(i));
end
y.an = aJ ;
y.bn = bJ ;
y.SigmaQ = sqrt(sigmaJ2);

