function [Yt,eht,Zt]=FIEGARCH_core(Ydata,d,omega,aL,bL,theta,delta, ErrorType,zeta);
%disp('[Yt eht Zt]=FIEGARCH_core(Ydata,d,omega,aL,bL,theta,delta, ErrorType, zeta')

%%%%%%%%%%%%%
% IMPORTANT %
%%%%%%%%%%%%%
%Set truncation lag and number of discarded values in FIXED PARAMETERS part

%%%%%%%%%%%%%%%%
%   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 %
%%%%%%%%%%%%
% Generates series of conditional variances eht for series of realizations
% Yt assumed to be generated by FIEGARCH(q,d,p) process. The function
% returns Nx3 matrix [Yt eht Zt] of the original process, conditional
% variances and original process standardized by conditional variances

%p=2
%q=1

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   INPUT CHECK %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   If ErrorType is not specified, muz is set to default value, i.e.sample
%   mean. 
%   7 and 9 entries option:
switch nargin
    case 7
    ErrorType='sample';
    %disp('muz=sample mean')
    case 9
    otherwise 
    error('7 or 9 inputs only: [Yt eht Zt]=FIEGARCH_core(Ydata,Z0,eht0,omega,aL,bL,theta,delta, ErrorType,zeta)');
end
%%%%%%%%%%%%%
% ErrorType %
%%%%%%%%%%%%%
%   ErrorType= Error distribution
%   zeta= parameters in error distribution function

%   VALID ALTERNATIVES:
%       'sample'(7 parameters only)
%       'stdnormal' by default 
%       '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:  DATA                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Ydata= data assumed to be generated by FIEGARCH(q,d,p)*Nx1


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   INPUT:  MODEL PARAMETERS        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% [d,omega,aL,bL,theta,delta]= parameters in the FIEGARCH(q,d,p)model
% for fiegarch(1,d,2):lag polynomials are of first and second order
% respectively with:
% a1=1; not used as a variable
a2=aL(1);
b1=bL(1); 

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


% muz= mean of |Zt|, by default calculated as mean of the observations,
% when ErrorType is set, muz is calculated based on assumed Zt
% distribution:  

if strcmp(ErrorType,'sample')
        muz=mean(abs(Ydata));
        %disp('Zt N(0,1) distributed')

elseif strcmp(ErrorType,'stdnormal')
        muz=0.7979;
        %disp('Zt N(0,1) distributed')

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

    
%%%%%%%%%%%%%%%%%%%%%%%
%   FIXED PARAMETERS  %
%%%%%%%%%%%%%%%%%%%%%%%
%   N= length of Ydata
%   JH: truncation of ht lag polynomial (AR terms)
%   J: truncation lag for aproximation FIEGARCH(0,d,J), i.e # of lags of
%   g(t) (MA terms)
%   Tstar= number of discarded values (size of generated eht<N)
N=length(Ydata);

Tstar=0;% To be set
JH=100;% 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)'])

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                     Nothing else to be specified                       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   INITIALIZING (PRESAMPLE) VALUES   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Yt0=JH initializing values of Yt: 
%           Yt(t)=0 for  t<=0 
Yt0=zeros(JH,1);
%   eht0=JH initializing values of ht:
%           eht(t)=sample variance of Y(t);  t=-JH+1,...,0   *JHx1
eht0=repmat(var(Ydata),JH,1);

Yt=[Yt0; Ydata];        % JH+Nx1
eht=[eht0; zeros(N,1)];    % JH+Nx1 
T=size(Yt,1);               % T=JH+N                   

pivec=cumprod(((0:1:J-3)-d)./(1:1:J-2));  % 1xJH-1: (1-L)^d coefficients
ARFIcoef=b1*[1 pivec(1:J-3)]-pivec; % (1-b1*L)(1-L)^d = 1-B1*L-...-BN*L^N



for t=JH+1:T;
  if eht(t-1)>0 && eht(t-1)<inf;
     
    Z12=Yt(t-(1:2))./sqrt(eht(t-(1:2))); % 2x1: [z(t-1),z(t-2)]'
  
    MA=[1 a2]*( theta*Z12 + delta*( abs(Z12)-muz ) );%(1x2)x(2x1)=1x1
  
    ARF=ARFIcoef*(log(eht(t-(1:JH)))-omega*ones(JH,1));  % (1xJH)x(JHx1)=(1,1) ARFIcoef*past logvariances

    eht(t)=exp(omega+ARF+MA); %variance
  else
      eht(t)=eht0(1);
      eht(t-1)=eht0(1);
       
      % break 
       
  end
 
end;

% Truncation if Tstar>0 chosen
if length(eht)==JH+N
    if N>Tstar
     %display ([num2str(Tstar), ' values discarded, (see fixed parameters: Tstar)'])
     Yt=Yt(JH+1+Tstar:end);          % N-Tstar x1 original data
     eht=eht(JH+1+Tstar:end);        % N-Tstar x1 generated conditional variance
     Zt=Yt(N-(0:N-1-Tstar))./sqrt(eht(N-(0:N-1-Tstar))); % Nx1 Yt standardized by eht
    else
     %display ('Number of discarded values > size of Ydata, no truncation applied (Tstar=0)')
     Yt=Yt(JH+1:end);         % Nx1 original data
     eht=eht(JH+1:end);       % Nx1 generated conditional variance
     Zt=Yt(N-(0:N-1))./sqrt(eht(N-(0:N-1))); % Nx1 Yt standardized by eht
    end
else
    %display ([num2str(Tstar), ' values discarded, (see fixed parameters: Tstar)'])
    %display (' NaN in generated volatility')
    Yt=0;           % default error value
    Zt=0;           % default error value
end
end



