%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   FIEGARCH SIMULATION    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [Y,Ht,Phi,Jumps]=FIEGARCH_Jump_simulate(N,d,omega,aL,bL,theta,delta,ErrorType,zeta)
%disp('**[Y,ht,Phi]=FIEGARCH_simulate(N,d,omega,aL,bL,theta,delta,ErrorType,zeta)')
%%%%%%%%%%%%%
% IMPORTANT %
%%%%%%%%%%%%%
% Set truncation lag and number of discarded values in FIXED PARAMETERS
% part lines 134-135
% Set Jump regime: Line 220

%%%%%%%%%%%%%%%%
%   COPYRIGHT  %
%%%%%%%%%%%%%%%%
%  Based on code for paper: Perez, Zaffaroni(2008),Finite-sample Properties of Maximum Likelihood and
%   Whittle Estimators in EGARCH and FIEGARCH Models, received on request from Ana Perez
%   25/01/2013

%   Written by:
%   Lucie Kraicova
%   kraicova.l@seznam.cz

%%%%%%%%%%%%
%   USAGE %
%%%%%%%%%%%%
% Simulate N observations from FIEGARCH(q,d,p) process using infinite
% FIEGARCH(0,d,infty) representation(eq. (11) Bollerslev & Mikkelsen
% (1996))truncated at lag J, effectively FIEGARCH(0,d,J).

%p=2
%q=1

%%%%%%%%%%%%
%   MODEL  %
%%%%%%%%%%%%

%   Y=z(t)exp[0.5*h(t)]
%   h(t)=w+Phi(L)*g(t-1)
%   gZ=theta*Zt+delta*(abs(Zt)-muz)
     
%       w= a0*(1-L)^(-d)*[a(L)/b(L)]
%       Phi(L)=(1-L)^(-d)*[a(L)/b(L)]
%       (1-L)^d= 1-d*sum_{k=1}^{infty}( Gamma(k-d)*Gamma(1-d)^{-1}*
%                   *Gamma(k+1)^{-1}*L^{k} ) ; Gamma(x)= gamma function
%       a(L)= sum a_{i}*L^{i} ; i=1,...p; a_{1}=1
%       b(L)= 1-sum b_{i}*L^{i} ; i=1,...q 

%   RESTRICTIONS:
%       -0.5<d<0.5  
%       -1<a1<1
%       -1<b2<1             
%       theta finite 
%       delta finite
%       delta*theta ~= 0


%%%%%%%%%%%%%
% ErrorType %
%%%%%%%%%%%%%
%   ErrorType= Error distribution
%   zeta= parameters in error distribution function

%   VALID ALTERNATIVES:
%       'stdnormal' by default (7 parameters only)
%       'stdged'->zeta=[nu],            1=<nu nu=2 equiv. to stdnormal
%       'stdstudent'->zeta=[nu],        2<nu
%       'skewstudent'->zeta=[nu,lambda] 2<nu -1<lambda<1

%   For stdged:     zeta=nu
%                   nu=1 -> Laplace distribution
%                   nu=2 -> N(0,1)
%                   nu->infty -> U(0,1)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   INPUT CHECK AND DEFAULT VALUES FOR OPTIONAL PARAMETERS %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%   Entries options, zeta decomposition and muz calculation:

switch nargin
    case 7
        ErrorType='stdnormal';
        muz=0.7979;
        %disp('Zt N(0,1) distributed')
    case 9
         if strcmp(ErrorType,'stdnormal')
                 zeta=[];
                 muz=0.7979;
         elseif strcmp(ErrorType,'stdged')
                 nu=zeta(1);
                 if nu<1
                     error('invalid nu')
                 end
                 muz=gamma(2/nu)/sqrt(gamma(3/nu)*gamma(1/nu));%from Anna Perez Toolbox
         elseif strcmp(ErrorType,'stdstudent')
                 nu=zeta(1);
                 if nu<2 || nu==2
                     error('invalid nu')
                 end
                 muz=1/sqrt(pi/2); %approximation by Kevin Sheppard in MFL Toolbox
         elseif strcmp(ErrorType,'skewstudent')
                 nu=zeta(1);
                 if nu<2 || nu==2
                     error('invalid nu')
                 end
                 skew=zeta(2);
                 if abs(skew)>0.99
                     error('invalid lambda')
                 end
                 muz=1/sqrt(pi/2); %approximation by Kevin Sheppard in MFL Toolbox
         else disp('invalid ErrorType')
         end
    otherwise
        error('7 or 9 inputs only:[Y,ht,Phi]=FIEGARCH_simulate(N,d,omega,aL,bL,theta,delta,ErrorType,zeta), [ErrorType,zeta] optional, stdnormal by default');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   INPUT:  PARAMETERS IN THE FUNCTION   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%   [d,omega,aL,bL]in ht equation
%   [theta,delta,muz]: in gt equation, muz calculated based on ErrorType
%   N: number of realizations

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   FIXED SIMULATION PARAMETERS  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%   JH: truncation of ht lag polynomial (AR terms)
%   J: truncation lag for aproximation FIEGARCH(0,d,J)
%   Tstar: values to be discarded to avoid starting-up problems

Tstar=1000;% To be set
JH=1000;% To be set

J=JH+length(aL)+1;
%disp(['Truncated at lag ', num2str(JH), ' of h(t)'])
%disp(['Eqivalently at lag ', num2str(J), ' of g(t)'])
%disp(['Discarded first ', num2str(Tstar), ' simulated values, to avoid starting-up problems'])

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%              Nothing else to be specified except for jumps (220)       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   INITIALIZING VALUES   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%

% phi_vec: initializing coefficients of Phi polynomial in FIEGARCH(0,d,J)
% ht: initializing variable_logvariance 
% Y: initializing variable_Y(t)=z(t)exp[0.5 h(t)]

phi_vec=zeros(1,J);                    
ht=ones(Tstar+N,1);
Y=ones(Tstar+N,1);  


%%%%%%%%%%%%%%%%%
%   SIMULATION  %
%%%%%%%%%%%%%%%%%

%   Simulate J+T*+N observations z(t) (J+T*+N column vector) to enable J
%   lags for simulation of h(1),next T*+N-1 z(t)'s for simulation of
%   (h(2)-h(T*+N))(moving window)and 1 z(t) for simulation of last
%   observation of Y (Y(T*+N)=z(T*+N)exp[0.5 h(T*+N)]).

Zt_length=J+Tstar+N;
if strcmp(ErrorType,'stdged')
    Zt=gedrnd(nu,[Zt_length,1]);
elseif strcmp(ErrorType,'stdnormal')
    Zt=random('Normal',0,1,[Zt_length,1]);
elseif strcmp(ErrorType,'stdstudent')
    Zt=stdtrnd(nu,Zt_length,1);
elseif strcmp(ErrorType,'skewstudent')
    Zt=skewtrnd(nu,skew,Zt_length,1);
else
    error('Invalid ErrorType input')
end

% Compute J+T*+N column vector
% of{g(z(1-J)),...,g(z(0)),g(z(1)),...,g(z(T*+N)}
gZ=theta*Zt+delta*(abs(Zt)-muz);

% Define vector of Phi polynomial coefficients:

% Compute 1xJ-2 and 1xJ-1 vectors of coefficients {pi(1),...,pi(J-2)} and
% D{D(1)...D(J-1)} representing weights stemming from 
% multiplication of EGARCH(0,d,J) lag polynomial (a(L)/b(L))by (1-L)^(-d):

pi_vec=cumprod(((0:1:J-3)+d)./(1:1:J-2)); % 1xJ-2 vector
D=[((0:1:J-2)+d)./(1:1:J-1)];% 1xJ-1 vector

%disp('coeff.Pi(1:10)')
%pi_vec(1:10)

%disp('coeff.D(1:10)')
%D(1:10)

%   Define 1xJ vector Phi=[phi(1),phi(2),...phi(J)]
%   For FIEGARCH(1,d,2):
% a1=1; not used as a variable
a2=aL;
b1=bL;

phi_vec(1)=1;                 
phi_vec(2)=b1*phi_vec(1)+(D(1)+a2); 
for j=3:J;                     
    phi_vec(j)=b1*phi_vec(j-1)+(D(j-1)+a2)*pi_vec(j-2);
end;
phi_vec=[phi_vec(1:J)];

%   Compute h(t)

for i=1:Tstar+N;  
   ht(i)=omega+phi_vec(end:-1:1)*gZ(i:J-1+i);
end

% @@@@@@@ Simulate jumps:@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

% 1) Additive jumps not entering the conditional volatility equation
     %   Compute (T*+N)x1 column vector of Y(t)without Jumps
         Y=Zt(J+1:length(Zt)).*exp(0.5*ht);
      
      
      %1.1) Sum of intraday Jumps: e.g. 24*(60/5)= 288 observations (24h trading, 5min data)
          
          % Settings:
            
           % Set 1/2: Number of jumps: Poisson(poiss_lambda) distribution
                 poiss_lambda=0.007; %0.014 Mancini(2012)+Facebook,0.028 simulated(1),0.007 for extreme case of rare huge jumps
                                  
               
           % Set 2/2: Size distribution: Normal(jumps_mu,jumps_sigma)
                 jumps_mu=0; %
               
                %jumps_sigma=abs(mean(Y*10)); % Link Jump variance to mean return
                 jumps_sigma=3; % 0.25 (Mancini 2012), 0.1 Facebook, 0.2 Simulated (1),3 for extreme case of rare huge jumps
                                  
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        Days_num=length(Y);
        Jumps=zeros(Days_num,1);
        for i=1:Days_num  
            jumps_num=poissrnd(poiss_lambda);
            jumps_intraday=normrnd(jumps_mu,jumps_sigma,jumps_num,1); 
            Jumps(i)=sum(jumps_intraday);
        end

        %   Compute (T*+N)x1 column vector of Y(t)
        Y=Y+Jumps;
   
  
 %@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@            
             
%   Discard first T* values ==> Nx1 column vector
Y=Y(Tstar+1:Tstar+N);
Jumps=Jumps(Tstar+1:Tstar+N);
Ht=[ht(Tstar+1:Tstar+N)];
Phi=(phi_vec)';
end









