function val = forecast_sdtsm_Quarter(Y, Phi, omega, sige2, gamma, Xbar, sigtau2, lenJ, aJ, bJ, sigmasJ, rlb, h)
% Forecast ten-year forward rate for each quarter
% Forecast starts at 197609 and ends at 201509
% INPUT:
%       Y:              ten-year foreward rate  Tx1 
%       Phi:            3x3
%       aJ;sigmasJ:     scalars
%       bJ;             1x3
%       h:              predicting horizons, scalar
% 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(T,1);
resid = nan(T,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 = 60: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;
        if j == h-2
            musJ = aJ(12) + bJ(12,:)*EX(2:4);  % (K-1)x1
            z1_temp = (musJ-rlb)./sigmasJ(12);
            yhat1 = rlb + (musJ-rlb).*normcdf(z1_temp) + sigmasJ(12).*normpdf(z1_temp);
        end
        if j == h-1
            musJ = aJ(12) + bJ(12,:)*EX(2:4);  % (K-1)x1
            z1_temp = (musJ-rlb)./sigmasJ(12);
            yhat2 = rlb + (musJ-rlb).*normcdf(z1_temp) + sigmasJ(12).*normpdf(z1_temp);
        end
    end
    musJ = aJ(12) + bJ(12,:)*EX(2:4);  % (K-1)x1
    z1_temp = (musJ-rlb)./sigmasJ(12);
    yhat3 = rlb + (musJ-rlb).*normcdf(z1_temp) + sigmasJ(12).*normpdf(z1_temp);
    yhat(ii+h) = (yhat1+yhat2+yhat3)/3;
    resid(ii+h) = mean(Y(ii+h-2:ii+h)) - yhat(ii+h) ;
end

yhat(isnan(yhat)) = [];
resid(isnan(resid)) = [];


rmsfe = sqrt(sum(resid.^2)/length(resid));

val.rmsfe = rmsfe ;
val.resid = resid;
val.yhat = yhat;


end