function val = forecast_jsz3(FR, yields, mats, h)
% Forecast ten-year foreward rate
% INPUT:
%       FR:             forward rates        Tx1  10y
%       yields:         yields               Txm  3m,6m,1-10y
%       mats:           maturities           mx1  3m,6m 12m-120m
%       h:              scalar predicting horizons
% OUTPUT:
%       RMSE
%       resid:  (T-h)x1  residuals from forecast
%       yhat:   (T-h)x1  fitted y


T = size(yields,1);
n = T - 40;        % total number of forecasts
yhat = nan(n,1);
resid = nan(n,1);
yields = yields/1200;
% AcP = AcP*1200;
% BcP = BcP*1200;

% Yin = yields(1:500,:);
% [eigv, ~] = eig(cov(Yin));
% eigv = fliplr(eigv) ;
% W = transpose(eigv) ;        % rows contain loadings, first row is first PC
% W(1,:) = W(1,:) / sum(W(1,:));
% W = W*1200  ;
% WN = W(1:3,:);

for ii = 20:n
    if mod(ii,50) == 0
        disp(ii);
    end
    Yin = yields(1:ii,:);
    [eigv, ~] = eig(cov(Yin));
    eigv = fliplr(eigv) ;
    W = transpose(eigv) ;        % rows contain loadings, first row is first PC
    W(1,:) = W(1,:) / sum(W(1,:));
    W = W*1200  ;
    EX = W(1:3,:)*yields(ii,:)';
    jszpars = estJSZ(Yin, W, 3, mats,[]) ;
    for j = 1:h
        EX = jszpars.KP_0P + jszpars.KP_PP * EX;
    end
    GATSM_loadings = GATSM(jszpars.KQ_0P, jszpars.KQ_PP, jszpars.Omega_cP, jszpars.rho0_cP, jszpars.rho1_cP, mats, 1);
    yhat(ii) = 1200*(GATSM_loadings.an(end) + GATSM_loadings.bn(end,:) * EX);
%     resid(ii) = yields(ii+h, 12) - yhat(ii);
    resid(ii) = FR(ii+h) - yhat(ii);
end

% resid = resid*1200;
rmsfe = sqrt(sum(resid(20:end).^2)/(n-19)) ;

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


end