%% Generate date from BPP model %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function data_dgp = dgp_FD(dgp_prmtr1,N,T,shocks_v,shocks_u,shocks_eps,shocks_eta)

%%
theta_dgp = dgp_prmtr1(1); 
gamma_eta_dgp = dgp_prmtr1(2); 
gamma_eps_dgp = dgp_prmtr1(3); 

N_dgp = N; 
T_dgp = T;
%%
data_dgp = zeros(N_dgp*T_dgp,2);
eps_t_dgp = zeros(N_dgp*T_dgp,1);
tao_dgp = zeros(N_dgp*T_dgp,1);
kappa_dgp = zeros(N_dgp*T_dgp,1);
%% 
 i_dgp = 1;
 
    while i_dgp <= N_dgp
    
    j_dgp = (i_dgp-1)*T_dgp+1;
    k_dgp = 1;
    
        while j_dgp <= i_dgp*T_dgp    
            
         if (k_dgp==1)
             
         shocks_eps_t = selif(shocks_eps(2,1:N)',shocks_eps(2,1:N)' ~= 999);
         eps_dgp = shocks_eps_t(ceil(size(shocks_eps_t,1)*rand(1,1)),1);
         lagged_eps_dgp = shocks_eps_t(ceil(size(shocks_eps_t,1)*rand(1,1)),1);
           
         shocks_eta_t = selif(shocks_eta(2,1:N)',shocks_eta(2,1:N)' ~= 999);
         eta_dgp = shocks_eta_t(ceil(size(shocks_eta_t,1)*rand(1,1)),1);             

         shocks_v_t = selif(shocks_v(2,1:N)',shocks_v(2,1:N)' ~= 999);
         v_dgp = shocks_v_t(ceil(size(shocks_v_t,1)*rand(1,1)),1);
 
         shocks_u_t = selif(shocks_u(2,1:N)',shocks_u(2,1:N)' ~= 999);
         u_dgp = shocks_u_t(ceil(size(shocks_u_t,1)*rand(1,1)),1);         
             
         else
             
         shocks_eps_t = selif(shocks_eps(k_dgp,1:N)',shocks_eps(k_dgp,1:N)' ~= 999);
         eps_dgp = shocks_eps_t(ceil(size(shocks_eps_t,1)*rand(1,1)),1);
           
         shocks_eta_t = selif(shocks_eta(k_dgp,1:N)',shocks_eta(k_dgp,1:N)' ~= 999);
         eta_dgp = shocks_eta_t(ceil(size(shocks_eta_t,1)*rand(1,1)),1);
      
            if (k_dgp==10 || k_dgp==11 || k_dgp==12)
                shocks_v_t = selif(shocks_v(9,1:N)',shocks_v(9,1:N)' ~= 999);
                v_dgp = shocks_v_t(ceil(size(shocks_v_t,1)*rand(1,1)),1);
      
                shocks_u_t = selif(shocks_u(9,1:N)',shocks_u(9,1:N)' ~= 999);
                u_dgp = shocks_u_t(ceil(size(shocks_u_t,1)*rand(1,1)),1);
 
            else     
                shocks_v_t = selif(shocks_v(k_dgp,1:N)',shocks_v(k_dgp,1:N)' ~= 999);
                v_dgp = shocks_v_t(ceil(size(shocks_v_t,1)*rand(1,1)),1);
 
                shocks_u_t = selif(shocks_u(k_dgp,1:N)',shocks_u(k_dgp,1:N)' ~= 999);
                u_dgp = shocks_u_t(ceil(size(shocks_u_t,1)*rand(1,1)),1);
            end
         end
        
        eps_t_dgp(j_dgp,1) = eps_dgp;

            if (j_dgp==(i_dgp-1)*T_dgp+1)
                tao_dgp(j_dgp,1) = eta_dgp;
                kappa_dgp(j_dgp,1) = gamma_eps_dgp*eps_dgp+u_dgp;
                data_dgp(j_dgp,1) = tao_dgp(j_dgp,1)+eps_dgp+theta_dgp*lagged_eps_dgp;       
            elseif (j_dgp>(i_dgp-1)*T_dgp+1)
                tao_dgp(j_dgp,1) = tao_dgp(j_dgp-1,1)+eta_dgp;
                kappa_dgp(j_dgp,1) = kappa_dgp(j_dgp-1,1)+gamma_eps_dgp*eps_dgp+u_dgp; 
                data_dgp(j_dgp,1) = tao_dgp(j_dgp,1)+eps_dgp+theta_dgp*eps_t_dgp(j_dgp-1,1);
            end

        data_dgp(j_dgp,2) = kappa_dgp(j_dgp,1)+gamma_eta_dgp*tao_dgp(j_dgp,1)+v_dgp;
       
                j_dgp = j_dgp+1;
                k_dgp = k_dgp+1;
        
        end
        i_dgp = i_dgp+1;
     
    end
    
end