function y = Grp(X, tau, kinfQ, lamQ, Xbar, gamma, Phi, mats, dt)
% Compute risk premium
% \Sigma*\lambda_t = \lambda_0 + \lambda_1 * [tau_t; X_t']'
% \lambda_0 = (I_n - Phi)*Xbar  - mu_Q
% \lambda-1 = [(I_n - Phi)*gamma, Phi - PhiQ]

% T = size(X,1);
Z = [tau X]; % Tx4
lambda0 = (eye(3)-Phi)*Xbar - [kinfQ;0;0];  % 3X1
lambda1 = [(eye(3)-Phi)*gamma, Phi-diag(lamQ)]; % 3x4
rp = lambda0' + Z*lambda1';  % Tx3

% loadsings Bm for GATSM
M = length(mats) ; 
% Atemp = 0 ;
Btemp = zeros(3,1) ;
% Ay = nan(1,M) ;
By = nan(3,M) ;
curr_mat = 1 ;
% rho0d = 0;
rho1d = ones(3,1) ;
% K0dp = [kinfQ;0;0]';
K1dp = diag(lamQ);
for i = 1:mats(M)
%     Atemp = Atemp + K0dp * Btemp + 0.5*transpose(Btemp) * omega * Btemp - rho0d ;
    Btemp = K1dp * Btemp - rho1d ;
    
    if i == mats(curr_mat)
%         Ay(1, curr_mat) = -Atemp / maturities(curr_mat) ;
        By(:,curr_mat) =  Btemp / mats(curr_mat) ;
        curr_mat = curr_mat +1 ;
    end
end
% A = Ay/dt ;   
B = By/dt ;  % 3*m
rp2 = rp * B;  % Txm

y.rp = rp;
y.rp2 = rp2;


