function val = forecast_sdtsm_rt(Y, theta0, prior, mats, lenJ, dt, rlb, h)
% real time oos for sdtsm

T = size(Y,1);
n = T -40;        % total number of forecasts
yhat = zeros(n,length(h));
resid = zeros(n,length(h));

iii = 1;
hmax = max(h);

rng('default')
rng(616)

options0  =  optimset('fminunc');
options0  =  optimset(options0 , 'Display', 'off');
options0  =  optimset(options0 , 'Diagnostics', 'off');

for ii = 20:n
    if mod(ii,50) == 0
        disp(ii);
    end
    Yin = Y(1:ii,:);
    [theta1,fval] = fminunc(@(theta)llk_sdtsm(theta, Yin', prior, rlb, mats, dt),theta0,options0);
    % theta1 = theta0;
    pars1 = theta2pars(theta1,mats,dt);
    mu_Z = [0; (eye(3)-pars1.Phi)*pars1.Xbar] ;
    Phi_Z = [1 zeros(1,3); (eye(3)-pars1.Phi)*pars1.gamma pars1.Phi] ;
    EX = esefit(Yin, pars1.Phi, pars1.omega, pars1.sige2, pars1.gamma, pars1.Xbar, pars1.sigtau2, pars1.aJ, pars1.bJ, pars1.sigmasJ, rlb, lenJ);
    for j = 1:hmax
        EX = mu_Z + Phi_Z*EX;
        if j == h(iii)
            musJ = pars1.aJ(12) + pars1.bJ(12,:)*EX(2:4);  % (K-1)x1
            z1_temp = (musJ-rlb)./pars1.sigmasJ(12);
            yhat(ii,iii) = rlb + (musJ-rlb).*normcdf(z1_temp) + pars1.sigmasJ(12).*normpdf(z1_temp);
            resid(ii,iii) = Y(ii+h(iii),12) - yhat(ii,iii);
            iii = iii+1;
        end
    end
    iii = 1;
    % musJ = pars1.aJ(12) + pars1.bJ(12,:)*EX(2:4);  % (K-1)x1
    % z1_temp = (musJ-rlb)./pars1.sigmasJ(12);
    % yhat(ii) = rlb + (musJ-rlb).*normcdf(z1_temp) + pars1.sigmasJ(12).*normpdf(z1_temp);
    % resid(ii) = Y(ii+h,12) - yhat(ii);
end

rmsfe = sqrt(sum(resid(20:end,:).^2,1)/(n-19)) ;

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

end

