function output = P_CopulaForecast(current_tau, params,data,copula_type,n_step_ahead)
% currect_tau: scalar, 当前tau值；
% params:3 by 1, 已估计的参数；
% data: [u v] 
% copula_type: 1 = SJC; 2 = Normal;
% n_step_ahead: 向前预测的步数；

%output：预测值；
if copula_type == 1 %SJC copula
    w_U = params(1);beta_U = params(2); alpha_U = params(3);
    w_L = params(4);beta_L = params(5); alpha_L = params(6);
    output(1,1) = current_tau(1);
    output(1,2) = current_tau(2);
    for n = 2:n_step_ahead+1
        if n <= 12
            x_U = w_U + beta_U*output(n-1,1) + alpha_U*mean(abs(data(end-12+n,1)-data(end-12+n,2)));
            x_L = w_L + beta_L*output(n-1,2) + alpha_L*mean(abs(data(end-12+n,1)-data(end-12+n,2)));
        else
            x_U = w_U + beta_U*output(n-1,1) ;
            x_L = w_L + beta_L*output(n-1,2) ;
        end
            output(n,1) = 0.998/(1+exp(x_U)) + 0.001;
            output(n,2) = 0.998/(1+exp(x_L)) + 0.001;
    end 
   output = output(2:end,:);
end

if copula_type == 2 %Gaussian copula
  % x = norminv(data(:,1),0,1);
  % y = norminv(data(:,2),0,1);
  x=data(:,1);
  y=data(:,2);
   w = params(1);beta = params(2); alpha = params(3);
   output(1) = current_tau; 
   for n = 2:n_step_ahead+1
       if n <=12
           psi = w + beta*output(n-1) + alpha*mean(x(end-12+n:end).*y(end-12+n:end));
       else
           psi = w + beta*output(n-1);
       end
           output(n) = 1.998/(1+exp(-psi))-0.999;
   end
   output = output(2:end)';
end


