## Generate date from BPP model ####################################################

dgp_FD <-  function(prmtr){

theta_dgp = prmtr[1]
gamma_eta_dgp = prmtr[2]
gamma_eps_dgp = prmtr[3]

N_dgp = N
T_dgp = Ts         

data_dgp = matrix(0,N_dgp*T_dgp,2)
eps_t_dgp = matrix(0,N_dgp*T_dgp,1)
tao_dgp = matrix(0,N_dgp*T_dgp,1)
kappa_dgp = matrix(0,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 <- t(shocks_eps[2, shocks_eps[2,] != 999])
         eps_dgp <- shocks_eps_t[1,ceiling(ncol(shocks_eps_t)*runif(1))]
         lagged_eps_dgp <- shocks_eps_t[1,ceiling(ncol(shocks_eps_t)*runif(1))]
         
         shocks_eta_t <- t(shocks_eta[2, shocks_eta[2,1:N] != 999])
         eta_dgp <- shocks_eta_t[1,ceiling(ncol(shocks_eta_t)*runif(1))]
 
         shocks_v_t = t(shocks_v[2, shocks_v[2,] != 999])
         v_dgp = shocks_v_t[1,ceiling(ncol(shocks_v_t)*runif(1))]
                
         shocks_u_t = t(shocks_u[2, shocks_u[2,] != 999])
         u_dgp = shocks_u_t[1,ceiling(ncol(shocks_u_t)*runif(1))]         
                 	       	
        }else{
        	
         shocks_eps_t <- t(shocks_eps[k_dgp, shocks_eps[k_dgp,] != 999])
         eps_dgp <- shocks_eps_t[1,ceiling(ncol(shocks_eps_t)*runif(1))]
         
         shocks_eta_t <- t(shocks_eta[k_dgp, shocks_eta[k_dgp,1:N] != 999])
         eta_dgp <- shocks_eta_t[1,ceiling(ncol(shocks_eta_t)*runif(1))]
         
         if (k_dgp ==10 || k_dgp==11 || k_dgp==12){
                shocks_v_t = t(shocks_v[9, shocks_v[9,] != 999])
                v_dgp = shocks_v_t[1,ceiling(ncol(shocks_v_t)*runif(1))]
                               
                shocks_u_t = t(shocks_u[9, shocks_u[9,] != 999])
                u_dgp = shocks_u_t[1,ceiling(ncol(shocks_u_t)*runif(1))]
                      
         }else { 
         	
         	shocks_v_t = t(shocks_v[k_dgp, shocks_v[k_dgp,] != 999])
            v_dgp = shocks_v_t[1,ceiling(ncol(shocks_v_t)*runif(1))]
                
            shocks_u_t = t(shocks_u[k_dgp, shocks_u[k_dgp,] != 999])
            u_dgp = shocks_u_t[1,ceiling(ncol(shocks_u_t)*runif(1))]
            
         }}

        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
            } else if (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]
            }
            
        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 of the inner loop
    i_dgp = i_dgp+1
    } # End of the outer loop
    
  return(data_dgp)    
}  # End of the dgp function


