% This program is to run the Monte Carlo simulations in Section 6-1.
% 
% There are two parameter sets to test alternatively: 
% 1. Fix r and run over delta as in Wn=n^{-delta}
% 2. Fix delta and run over r as in r^{-n}
% 
clear
rng shuffle
% --- Parameter Settings ---
TestOnDelta=0;              % set 1 to fix r, set others to fix delta
M=10000;                    % Monte Carlo runs
SZ=[100 200 300 400 500];   % sample size
h1=[1 1.2 0.32]; muu=1;
h2=[1 0.5 0.06]; muz=2;
h3=[1 0.9 0.8];             
h4=[1 0.7 0.12]; mup=2;
nn=200;                     % Riemann sum intervals for MISE
alpha_l=sum(h3);
r_chopoff=0.025;            % chop off rate on both sides, e.g. 5%
ell=3;                      % define the ell lags
% --- End of Parameter setting ---
if TestOnDelta==1
    delta=0.01:0.01:0.9;    % run over delta
    r=0.9;                  % fix r
    K=length(delta);
else
    delta=0.19;             % fix delta
    r=0.1:0.01:1;           % run over r
    K=length(r);
end
NN=length(SZ);              % how many sample sizes to test
MISE=zeros(NN,K);           % matrix to store the MISE results
mise=zeros(1,M);            % a temporary variable
% --- Monte Carlo simulations start ---
for k=1:NN                  % go through various sample sizes
    N=SZ(k);
for j=1:K                   % go through various delta(s) or r(s)
    if TestOnDelta==1       % on-line report during the runs
        fprintf('N=%d --- r=%4.2f --- delta=%4.2f\n',N,r,delta(j))
    else
        fprintf('N=%d --- r=%4.2f --- delta=%4.2f\n',N,r(j),delta)
    end
    
    for i=1:M               % M Monte Carlo runs
        % --- Generate I/O Data ---
        e=randn(N,1);
        Un=muu+filter(h1,1,e);
        Wn=sign(Un).*((abs(Un)+1).^2-1);
        deltan=unifrnd(-1,1,N,1);
        Zn=muz+filter(h2,1,deltan);
        Vn=Zn+Wn;
        Xn=filter(h3,1,Vn);
        deltan=unifrnd(-1,1,N,1);
        Pn=mup+filter(h4,1,deltan);
        Yn=Xn+Pn;
        % --- End of Generate I/O Data ---
        % --- Calculate MISE ---
        CI=[quantile(Un,r_chopoff) quantile(Un,1-r_chopoff)];% Confidence Interval
        u=linspace(CI(1),CI(2),nn); % nn equal length of points in CI
        if TestOnDelta==1           % run over delta or r
            ml_u_hat=semiProdGaussianX(u,Un,Yn,ell,delta(j),r);
        else
            ml_u_hat=semiProdGaussianX(u,Un,Yn,ell,delta,r(j));
        end
        mu=sign(u).*((abs(u)+1).^2-1);
        [tmp,I]=min(abs(mu));
        ml_0_hat=ml_u_hat(I);
        du=(u(2)-u(1));
        mise(i)=sum(((ml_u_hat-ml_0_hat)/alpha_l-mu).^2)*du;
    end
    MISE(k,j)=mean(mise);   % take average of M MISE(s)
end
end
% --- End of Monte Carlo simulations ---

% --- prepare to save the MISE results and make a preliminary plot
t=clock;
dirpath='';         % set directory of where to save the results;
if TestOnDelta==1   % delta or r
    plot(delta,MISE','LineWidth',2)
    xlabel('\delta')
    [minMISE, Idx]=min(MISE,[],2);
    text(delta(Idx),minMISE,'X')
    str=strcat('Estar_delta_r_',num2str(r*100),'_',num2str(t(2)),num2str(t(3)),num2str(t(4)),num2str(t(5)));
    str=fullfile(dirpath,str);
    save(str, 'MISE','delta','r','SZ')
else
    plot(r,MISE','LineWidth',2);
    xlabel('r')
    [minMISE, Idx]=min(MISE,[],2);
    text(r(Idx),minMISE,'X')
    str=strcat('Estar_r_delta_',num2str(delta*100),'_',num2str(t(2)),num2str(t(3)),num2str(t(4)),num2str(t(5)));
    str=fullfile(dirpath,str);
    save(str, 'MISE','delta','r','SZ')
end
title('E^*')
ylabel('MISE')
legend(string(SZ))
grid minor
% --- end of program ---

function m_hat=semiProdGaussianX(u,Un,Yn,ell,delta,r)
% The semirecursive and nonparametric algorithm (3.6)
% u:    new input data for prediction  
% Un:   measured inputs
% Yn:   measured outputs
% ell:  input lags to used
% delta:parameter that determines window bandwidth
% r:    parameter ralated to convergence rate
    if size(u,1)==1
        u=u';       % convert to a M x 1 vector if not
    end
    N=length(Un);   % size of measured data
    M=length(u);    % size of new data
    K=@(x) exp(-0.5*sum(x.^2,2));% product Gaussian kernel
    Wn=(1:N).^(-delta);

    if r==2         % E_2 in Section 5, a special case of E
        rn=Wn.^ell./cumsum(Wn.^ell);
    else            % r<1 for E and r=1 for E_1, a special case of E 
        rn=(1:N).^(-r);
    end
    gn=zeros(M,1); fn=gn;   % initial values for the semi-recursive algorithm
    for i=1:N               % go through measured data
        if (i-ell+1)<1
            s=i:-1:1;
        else
            s=i:-1:(i-ell+1);
        end
        TMP=(u*ones(1,length(s))-ones(M,1)*Un(s)')/Wn(i); % M x ell
        gn=(1-rn(i))*gn + rn(i)*Yn(i)*K(TMP)/Wn(i)^ell;
        fn=(1-rn(i))*fn + rn(i)*K(TMP)/Wn(i)^ell;
    end
    m_hat=gn./fn;
    if size(m_hat,2)==1
        m_hat=m_hat';
    end