function [INS] = Func_PLS_INS(Ret_ac, Ret_x, tdate, zsc1, cons1, ma1)
% tdate: longest time period (1985 - 2015)
% zsc1: z = normalization in 1st step, otherwise no zscore (in-sample)
% cons: 0 = no intercept in the 2nd step, 1=with intercept 
% ma: 0 = no moving average in the 3rd step, 12=12-month moving average

beta = []; tstat = []; R2ins = [];
INS = [];

for s = min(Ret_ac(:,6)):max(Ret_ac(:,6))
    Ret_s = Ret_ac(Ret_ac(:,6)==s,:);
    Y = Ret_s(:,2);
    tY = Ret_s(:,1);
    
    Pi_s = [];
    RetN = [];
   for i = min(Ret_x(:,6)):max(Ret_x(:,6))
        x_i  = nan*tdate;                                    %% Construct a big matrix
        Ret_i = Ret_x(Ret_x(:,6)==i,:);
        tX = Ret_i(:,1); 
        X = Ret_i(:,7);        
        if zsc1 == 'z'
            X_ins = zscore(X);
            Y_ins = Y-mean(Y);
        else
            X_ins = X;
            Y_ins = Y;
        end
           if tY(1) <= tX(1) && tY(end) >= tX(end)
            tt = find(tY(:)==tX(1));
            td = find(tY(:)==tX(end));
            rh = Y_ins(tt+1:td);                                     
            lh = X_ins(1:end-1);                                    
            Pi_s = [Pi_s, olsnw(lh,rh)];                                        
           
           elseif tY(1) <= tX(1) && tY(end) < tX(end)
                tt = find(tY(:)==tX(1));
                td = find(tX(:)==tY(end));
                rh = Y_ins(tt+1:end);
                lh = X_ins(1:td-1);
                Pi_s = [Pi_s, olsnw(lh,rh)];
 
            elseif tY(1)> tX(1) && tY(end) <= tX(end)
                tt = find(tX(:)==tY(1));
                td = find(tX(:)==tY(end));
                rh = Y_ins(1:end);
                lh = X_ins(tt-1:td-1);
                Pi_s = [Pi_s, olsnw(lh,rh)];

           elseif tY(1)> tX(1) && tY(end)> tX(end)
                tt = find(tX(:)==tY(1));
                td = find(tY(:)==tX(end));
                rh = Y_ins(1:td+1);
                lh = X_ins(tt-1:end);
                Pi_s = [Pi_s, olsnw(lh,rh)];
           end
                ta = find(tdate(:)==tX(1));
                tb = find(tdate(:)==tX(end));
                x_i(ta:tb) = X_ins(1:end);
                RetN = [RetN, x_i];         
   end
     ds = nan(length(tY),1);
   if cons1 == 0
       for t = 1:length(tY)
          RetN_t = RetN(t,:)';
          aa = find(~isnan(RetN_t));
          Rt = RetN_t(~isnan(RetN_t),:);                              
          Pi_sT = Pi_s(2,:)';                                              
          rhv = Pi_sT(aa);
          bv = olsnw(Rt, rhv, 0);
          ds(t,1) = bv(1);
       end
   else
       for t = 1:length(tY)
          RetN_t = RetN(t,:)';
          aa = find(~isnan(RetN_t));
          Rt = RetN_t(~isnan(RetN_t),:);                              
          Pi_sT = Pi_s(2,:)';                                              
          rhv = Pi_sT(aa);
          bv = olsnw(Rt, rhv);
          ds(t,1) = bv(2);
       end
   end
   if ma1 == 0
        lhv = Y(2:end)*100;                                                    
        rhs = [ones(size(ds(1:end-1))),  zscore(ds(1:end-1))];   
        [beta_s, tstat_s, s2_s, vcvnw_s, R2_s] = olsnw(lhv,rhs,0);
   else 
        lhv = Y(13:end)*100;                                                                       
        ma_ds = tsmovavg(ds, 'm', 12, 1);                
        rhs = [ones(size(ma_ds(12:end-1))),  zscore(ma_ds(12:end-1))]; 
        [beta_s, tstat_s, s2_s, vcvnw_s, R2_s] = olsnw(lhv,rhs,0);
   end
       beta = beta_s(2);
       tstat = tstat_s(2);
       R2ins = R2_s*100;
       INS = [INS; [beta, tstat, R2ins]];
end

end