function val = forecast_sdtsm2(Y, Phi, omega, sige2, gamma, Xbar, sigtau2, lenJ, aJ, bJ, sigmasJ, rlb, h)
% Forecast ten-year forward rate
% INPUT:
%       yields:         ten-year foreward rate  Tx1
%       Phi:            3x3
%       aJ;sigmasJ:     scalars
%       bJ;             1x3
%       h:              scalar predicting horizons
% OUTPUT:
%       RMSE
%       resid:  (T-h)x1  residuals from forecast
%       yhat:   (T-h)x1  fitted y


T = size(Y,1);
n = T -30;        % total number of forecasts
yhat = nan(n,1);
resid = nan(n,1);
mu_Z = [0; (eye(3)-Phi)*Xbar] ;
Phi_Z = [1 zeros(1,3); (eye(3)-Phi)*gamma Phi] ;

if any(abs(eig(Phi))>1)
    val = nan;
    return
end
Sigma_tilde = omega - (gamma * gamma')*sigtau2 ;
if invalidCovMat(Sigma_tilde)
    val = nan ;
    return
end


for ii = 58:n
%     EX = [X(ii,:), tau(ii)]';  % 4x1
    EX = esefit(Y(1:ii,:), Phi, omega, sige2, gamma, Xbar, sigtau2, aJ, bJ, sigmasJ, rlb, lenJ);
    for j = 1:h
        EX = mu_Z + Phi_Z*EX;
    end
    musJ = aJ(12) + bJ(12,:)*EX(2:4);  % (K-1)x1
    z1_temp = (musJ-rlb)./sigmasJ(12);
    yhat(ii) = rlb + (musJ-rlb).*normcdf(z1_temp) + sigmasJ(12).*normpdf(z1_temp);
    resid(ii) = Y(ii+h,12) - yhat(ii);
end

rmsfe = sqrt(sum(resid(60:end).^2)/(n-59));

val.rmsfe = rmsfe ;
val.resid = resid(58:end);
val.yhat = yhat;


end
