function FIEGARCH_montecarlo_final(scenario)
display('FIEGARCH_montecarlo_final(scenario)')

%! IMPORTANT:   Set truncation lags in FIEGARCH_simulate and FIEGARCH_core
%               Set parameter bounds LB,UB,LBMLE,UBMLE
%               Boundary conditions + dealing with boundary coefficients
%               (in WWE)+dealing with jumps in WWE and W and jump settings
%               in FIEGARCH_Jump_Simulate, choose FIEGARCH_simulate without
%               jumps line 564 and lock jump detectionin WWE and W
%               Number of estimations with measured time
%               Check if desired estimators unlocked (using '%{,%}'in the code)
%               Add a label to the results at the very bottom of the code

         
%%%%%%%%%%%%
%   USAGE  %
%%%%%%%%%%%%
% Perform MonteCarlo Experiments on the finite sample properties of the
% Maximum lilelihood, Whittle and Wavelet Whittle estimators of parameters
% of FIEGARCH(1,d,2) process, where parameters are subject to constrains.
% See Bollerslev & Mikkelsen (1996, 1999). Various distributions of errors allowed.
%The estimation uses simple minimization using fmincon without using user-supplied gradient, 
% optimization of the setup would improve the estimates.

%%%%%%%%%%%%%%%%
%   COPYRIGHT  %
%%%%%%%%%%%%%%%%
%   Lucie Kraicova
%   kraicova.l@seznam.cz

%%%%%%%%%%%%%%%%%%%%%%%%%
%   SCENARIOS           %
%%%%%%%%%%%%%%%%%%%%%%%%%

if  scenario==1
    M=1000;  %num of simulations
    
    N=2^11; %11 2048 % length of simul data
            %12 4096
            %14 16384 no MLE
            %16 65536 no MLE
    
    Errortype= 'default'; %if default, zeta is irrelevant
    %zeta=0 ;
    
    d=0.25; % set all param
    aL=0.5;
    bL=0.5;
    theta=-0.3;
    delta=0.5;
    
    method1='modwt';  % currently WWE1 only - unlock other two if needed (FWE unlocked, MLE locked)
    method2='modwt';
    method3='DWT';
    
    wavelet1='d4';
    wavelet2='d4';
    wavelet3='la8';
    
    level1=9; % J-3
    level2=9;
    level3=9;
    
elseif scenario==2
    M=1000;    %1000
    N=2^14; %11 2048 
            %12 4096
            %14 16384 
            %16 65536 
    
    Errortype= 'default'; %if default, zeta is irrelevant
    zeta=0 ;
    
    d=0.25;
    aL=0.5;
    bL=0.5;
    theta=-0.3;
    delta=0.5;
    
    method1='modwt';
    method2='modwt';
    method3='DWT';
    
    wavelet1='d4';
    wavelet2='d4';
    wavelet3='la8';
    
    level1=12;
    level2=12;
    level3=7;
    
elseif scenario==3
    M=1000;    %1000
    N=2^11; %11 2048
            %12 4096
            %14 16384 
            %16 65536 
    
    Errortype= 'default'; %if default, zeta is irrelevant
    %zeta=0 ;
    
    d=0.45;
    aL=0.5;
    bL=0.5;
    theta=-0.3;
    delta=0.5;
    
    method1='modwt';
    method2='modwt';
    method3='DWT';
    
    wavelet1='d4';
    wavelet2='d4';
    wavelet3='la8';
    
    level1=9;
    level2=9;
    level3=8;
    
elseif scenario==4
    M=1000;    %1000
    N=2^9;  %11 2048
            %12 4096
            %14 16384 
            %16 65536 
    
    Errortype= 'default'; %if default, zeta is irrelevant
    %zeta=0 ;
    
    d=0.45;
    aL=0.5;
    bL=0.5;
    theta=-0.3;
    delta=0.5;
    
    method1='DWT';
    method2='DWT';
    method3='DWT';
    
    wavelet1='Haar';
    wavelet2='d4';
    wavelet3='la8';
    
    level1=9;
    level2=7;
    level3=6;
    
elseif scenario==5
    M=1000;    %1000
    N=2^9; %11 2048
            %12 4096
            %14 16384 
            %16 65536 
    
    Errortype= 'default'; %if default, zeta is irrelevant
    %zeta=0 ;
    
    d=0.45;
    aL=0.5;
    bL=0.5;
    theta=-0.3;
    delta=0.5;
    
    method1='DWT';
    method2='DWT';
    method3='DWT';
    
    wavelet1='Haar';
    wavelet2='d4';
    wavelet3='LA8';
    
    level1='dwtcut4';
    level2='dwtcut4';
    level3='dwtcut4';
 
elseif scenario==6
    M=1000;    %1000
    N=2^11; %11 2048
            %12 4096
            %14 16384 
            %16 65536 
    
    Errortype= 'default'; %if default, zeta is irrelevant
    %zeta=0 ;
    
    d=0.25;
    aL=0.5;
    bL=0.5;
    theta=-0.3;
    delta=0.5;
    
    method1='DWT';
    method2='DWT';
    method3='DWT';
    
    wavelet1='Haar';
    wavelet2='d4';
    wavelet3='LA8';
    
    level1=4;
    level2=4;
    level3=4;   
    
else
            error('invalid scenario')
end

% Setting Standard Normal distribution as default
if strcmp(Errortype,'default')
        ErrorType='stdnormal';
        zeta=[];
        nu='NaN';
        skew='NaN';
end

% Extensions
zetaflag=0; %Choose 1 if distribution parameters should be estimated

%omega=0; % Choose any real number as a conditional log-variance mean
omega=-7; % -7 based on Facebook logreturns time series

% Time measuring: for how many rounds (<=M)
Measure=5;


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


%   Entries options, zeta decomposition :


if strcmp(ErrorType,'stdnormal')
                 zeta=[];
                 nu='NaN';
                 skew='NaN';
                 
elseif strcmp(ErrorType,'stdged')
                 nu=zeta(1);
                 if nu<1
                     error('invalid nu')
                 end
                 skew='NaN';
                 
elseif strcmp(ErrorType,'stdstudent')
                 nu=zeta(1);
                 if nu<2 || nu==2
                     error('invalid nu')
                 end
                 skew='NaN';
                 
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 (skew)')
                 end
                 
else
        error('invalid ErrorType')
         
end



%%%%%%%%%%%%%%%%%%%%
%  INPUT: METHOD   %
%%%%%%%%%%%%%%%%%%%%
% Methods Available:
%   'DWT' (uses wmtsa toolbox)
%       Fixed parameters to be specified in FIEGARCH_L_WW_estimate and
%       FIEGARCH_L_WW_estimate_B:
%           boundary   -- (optional) boundary conditions to use (string)
%                   Valid values: 'circular' or 'reflection'
%                   Default: 'reflection'
%                   Chosen: 'circular'
%           * opts       -- (optional) Additional function options
%              Optional input arguments are specified as name-value pairs:
%                   * RetainVJ -- Boolean flag indicating whether to
%                           scaling coefficients have been retained at all
%                           levels.
%                       Values:  1 = true,  all VJ retained,
%                                0 = false, VJ retained for J0 level.
%                       Default: 0, VJ retained only at J0 level.
        

%   'dwt' (uses matlab toolbox (function wavedec))


%%%%%%%%%%%%%%%%%%%%
%  INPUT: WAVELETS %
%%%%%%%%%%%%%%%%%%%%
% Filters Available:
%   For WMTSA transform(MODWT): write: 
%       wtfilter('list')
%       or a valid filter struct

%   For MATLAB transform(DWT,SWT,NDWT): write: 
%       waveletfamilies('n')
%       or a valid filter struct

%   pay attention when choosing filters: e.g. Haar~D2(WMTSA)~db1(MATLAB)

%%%%%%%%%%%%%%%%%%%%
%   INPUT: LEVEL   %
%%%%%%%%%%%%%%%%%%%%
%   0 < integer < log2(T)
%   'conservative' <log2( (T / (L - 1)) - 1) (L is the wavelet support)
%   'max' <log2(T)
%   'supermax' <log2(1.5 * T)

%'REFLECTION'
%   'dwtcut3'  <floor(log2(T))-3  16 J-level coefficients in DWT
%   'dwtcut4'  <floor(log2(T))-4  32 J-level coefficients in DWT
%   'dwtcut5'  <floor(log2(T))-5  64 J-level coefficients in DWT
%   'dwtcut6'  <floor(log2(T))-6 128 J-level coefficients in DWT
%   'dwtcut7'  <floor(log2(T))-7 256 J-level coefficients in DWT
%   'dwtcut8'  <floor(log2(T))-8 512 J-level coefficients in DWT

% 'CIRCULAR'
%   'dwtcut3'  <floor(log2(T))-3  8 J-level coefficients in DWT
%   'dwtcut4'  <floor(log2(T))-4  16 J-level coefficients in DWT
%   'dwtcut5'  <floor(log2(T))-5  32 J-level coefficients in DWT
%   'dwtcut6'  <floor(log2(T))-6  64 J-level coefficients in DWT
%   'dwtcut7'  <floor(log2(T))-7 128 J-level coefficients in DWT
%   'dwtcut8'  <floor(log2(T))-8 256 J-level coefficients in DWT


%%%%%%%%%%%%%%%%%%%%%%%
%    PARAMETERS       %
%%%%%%%%%%%%%%%%%%%%%%%
% parameters=[omega,aL,bL,theta,delta,zeta]

%   [omega,aL,bL,theta,delta]= parameters in the FIEGARCH(q,d,p)model 
%   for FIEGARCH(1,d,2):lag polynomials (bL,aL) are of first and second
%   order respectively with:

% a1=1; %not used as a variable
a2=aL(1);
b1=bL(1); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   INPUT:  OTHER PARAMETERS        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   N number FIEGARCH realizations
%   M number of experiments
%   zetaflag 1 if zeta should be estimated, 0 if zeta fixed at zeta0

%   ErrorType= distribution of Zt
%   zeta=parameters in Zt distribution function
%       VALID ALTERNATIVES:
%       'stdnormal' by default (8 parameters only)or 10 with zeta=[],
%       'stdged'->zeta=[nu],            1=<nu nu=2 equiv. to stdnormal
%                                             nu=1 equivalent to Laplace distribution
%                                             nu->inf equiv. to uniform distribution
%       '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) Standard Normal distribution
%                       nu->infty -> U(0,1) Uniform distribution
  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                     Check bounds for estimates below                   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%INITIALIZING  PARAMETERS %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%parameters=[d0,omega0,aL0,bL0,theta0,delta0]
%set to true values by default

%{
d0=d;
omega0=omega;
a20=a2;          
b10=b1;                
theta0=theta;        
delta0=delta; 
%}

% Test convergence using zeros as initializing

d0=0;
omega0=0;
a20=0;          
b10=0;                
theta0=-0.01;        
delta0=0; 


% Restrictions:
    %  -0.5<d<0.5
    %  omega=0 assumption
	%  -1<a1<1  
    %  -1<b1<1             
    % 0<theta_2
    % delta finite
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   INITIALIZING  VALUES   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
d_vec_MLE=zeros(M,1);
omega_vec_MLE=zeros(M,1);
a2_vec_MLE=zeros(M,1);
b1_vec_MLE=zeros(M,1);
theta_vec_MLE=zeros(M,1);
delta_vec_MLE=zeros(M,1);
nu_vec_MLE=zeros(M,1);
skew_vec_MLE=zeros(M,1);


d_vec_W=zeros(M,1);
omega_vec_W=zeros(M,1);
a2_vec_W=zeros(M,1);
b1_vec_W=zeros(M,1);
theta_vec_W=zeros(M,1);
delta_vec_W=zeros(M,1);
nu_vec_W=zeros(M,1);
skew_vec_W=zeros(M,1);



d_vec_WW1=zeros(M,1);
omega_vec_WW1=zeros(M,1);
a2_vec_WW1=zeros(M,1);
b1_vec_WW1=zeros(M,1);
theta_vec_WW1=zeros(M,1);
delta_vec_WW1=zeros(M,1);
nu_vec_WW1=zeros(M,1);
skew_vec_WW1=zeros(M,1);


d_vec_WW2=zeros(M,1);
omega_vec_WW2=zeros(M,1);
a2_vec_WW2=zeros(M,1);
b1_vec_WW2=zeros(M,1);
theta_vec_WW2=zeros(M,1);
delta_vec_WW2=zeros(M,1);
nu_vec_WW2=zeros(M,1);
skew_vec_WW2=zeros(M,1);


d_vec_WW3=zeros(M,1);
omega_vec_WW3=zeros(M,1);
a2_vec_WW3=zeros(M,1);
b1_vec_WW3=zeros(M,1);
theta_vec_WW3=zeros(M,1);
delta_vec_WW3=zeros(M,1);
nu_vec_WW3=zeros(M,1);
skew_vec_WW3=zeros(M,1);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  RESTRICTIONS ON ESTIMATES   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Restrictions:
    %  -0.5<d<0.5
    %  omega=0 assumption (W,WW)
	%  -1<a1<1  
    %  -1<b1<1             
    % 0<theta_2: Bounds for estimation: only positive or negative!
    % delta finite

% Although theoretically infinitely high/ low values for some parameters, somewhat more
% reasonable intervals chosen
LB=[-0.48 0  -0.9 -0.9 -0.9 -0.9]; % Lower bound (d,omega,a2,b1,theta,delta)
UB=[ 0.48 0   0.9  0.9  0  0.9]; % Upper bound (d,omega,a2,b1,theta,delta)
parameters0=[d0,omega0,a20,b10,theta0,delta0];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   DISTRIBUTION OF Zt: PARAMETERS AND RESTRICTIONS %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if zetaflag==1
    if strcmp(ErrorType,'stdnormal')
                %LB=LB;
                %UB=UB;           
                %zeta0=[];
                %parameters0=[parameters0,zeta0];
                
                
    elseif strcmp(ErrorType,'stdged')
                 
                 LB=[LB 1];    % Lower bounds for  (d,omega,a2,b1,theta,delta,nu)
                 UB=[UB 10];   % Upper bounds for  (d,omega,a2,b1,theta,delta,nu)                
                 zeta0=[nu];
                 parameters0=[parameters0,zeta0];
                 
                 
    elseif strcmp(ErrorType,'stdstudent')
                 LB=[LB 2.02];    % Lower bounds for  (d,omega,a2,b1,theta,delta,nu)
                 UB=[UB 10];    % Upper bounds for  (d,omega,a2,b1,theta,delta,nu)
                 zeta0=[nu];
                 parameters0=[parameters0,zeta0];
                 
                 
    elseif strcmp(ErrorType,'skewstudent')
                 LB=[LB 2.02 -0.98];    % Lower bounds for  (d,omega,a2,b1,theta,delta,nu,skew)
                 UB=[UB 10 0.98];     % Upper bounds for  (d,omega,a2,b1,theta,delta,nu,skew)
                 zeta0=[nu,skew];
                 parameters0=[parameters0,zeta0];
                
    else disp('invalid ErrorType')
    end
end

% Allow MLE to estimate omega
LBMLE=[ LB(1) -20 LB(3:length(LB))]; % Lower bound (d,omega,a2,b1,theta,delta(,nu,skew))
UBMLE=[ UB(1)  20 UB(3:length(UB))]; % Upper bound (d,omega,a2,b1,theta,delta(,nu,skew))

%disp(['ErrorType=', num2str(ErrorType)])
%disp(['parameters0=', num2str(parameters0)])
%disp(['LB=', num2str(LB)]);
%disp(['UB=', num2str(UB)]);

%options=optimset('LargeScale','off','Display','off','TolFun',1e-4); %perez
%MLE
options=optimset('LargeScale','off','Algorithm','active-set','Display','off');
%options=optimset('LargeScale','on','Display','off'); %active-set when
%gradient not provided

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    INITIALIZE DIARY RECORD    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

diary Diary_file %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
diary on;
disp(' Jumps '); % are jumps enabled
disp('Montecarlo: start');
t0=clock;
disp(datestr(clock,0));
disp('  ')
diary off;

%%%%%%%%%%%%%%%%%%%%%%%%%%
%    START MONTECARLO    %
%%%%%%%%%%%%%%%%%%%%%%%%%%


if zetaflag==0
    
    %for i=1:M ; 
    
    matlabpool open % lock or unlock both open here and close at the bottom  
    parfor i=1:M
        
        disp(['i=', int2str(i)]);
        
        % Simulate data
        % A) No Jumps
             %[Ydata,HtSim,PhiSim]=FIEGARCH_simulate(N,d,omega,aL,bL,theta,delta,ErrorType,zeta); % Nx1
        
        % B) Jumps
        % Adjust jump settings in the function
             [Ydata,HtSim,PhiSim,Jumps]=FIEGARCH_Jump_simulate(N,d,omega,aL,bL,theta,delta,ErrorType,zeta); % Nx1
        
        tic

    %{ 
        % MLE
        %Measure time
        if i<=Measure
         diary Diary_file %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
         diary on;
         disp(datestr(clock,0));
         disp('Start MLE');
         toc
         tic
         diary off;
        end
        
       
        % Minimize negative log-likelihood
        f_MLE=@(parameters)FIEGARCH_L_ML_estimate(parameters,Ydata,ErrorType,zeta);
        parametershat_MLE=fmincon(f_MLE,parameters0,[],[],[],[],LBMLE,UBMLE,[],options);
        
        % Store estimates in vectors
        d_vec_MLE(i)=parametershat_MLE(1);
        omega_vec_MLE(i)=parametershat_MLE(2);
        a2_vec_MLE(i)=parametershat_MLE(3);       
        b1_vec_MLE(i)=parametershat_MLE(4);              
        theta_vec_MLE(i)=parametershat_MLE(5);   
        delta_vec_MLE(i)=parametershat_MLE(6);  
       
        % Generate vectors for fixed parameters
        if strcmp(ErrorType,'stdged')
            nu_vec_MLE(i)=nu;
        elseif strcmp(ErrorType,'stdstudent')
            nu_vec_MLE(i)=nu;
        elseif strcmp(ErrorType,'skewstudent')
            nu_vec_MLE(i)=nu;
            skew_vec_MLE(i)=skew;
       
        end;
      
        
        %Measure time
        if i<=Measure
         diary Diary_file %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
         diary on;
         disp('End MLE, Start Whittle');
         disp(datestr(clock,0));
         toc
         tic
         diary off;
        end
        %}
       
        

    % Whittle
        % Minimize the Whittle objective function
        f_W=@(parameters)FIEGARCH_L_W_estimate(parameters,Ydata,ErrorType,zeta);
        parametershat_W=fmincon(f_W,parameters0,[],[],[],[],LB,UB,[],options);
        
        % Store estimates in vectors
        d_vec_W(i)=parametershat_W(1);
        omega_vec_W(i)=parametershat_W(2);
        a2_vec_W(i)=parametershat_W(3);       
        b1_vec_W(i)=parametershat_W(4);              
        theta_vec_W(i)=parametershat_W(5);   
        delta_vec_W(i)=parametershat_W(6);  
        
        % Generate vectors for fixed parameters
        if strcmp(ErrorType,'stdged')
            nu_vec_W(i)=nu;
        elseif strcmp(ErrorType,'stdstudent')
            nu_vec_W(i)=nu;
        elseif strcmp(ErrorType,'skewstudent')
            nu_vec_W(i)=nu;
            skew_vec_W(i)=skew;
        end;
       
        
        
        %Measure time
        if i<=Measure
         diary Diary_file %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
         diary on;
         disp(datestr(clock,0));
         disp('End Whittle, Start Wavelet Whittle');
         toc
         tic
         diary off;
        end

      
    % Wavelet Whittle 1
        % Minimize the Wavelet Whittle objective function
        f=@(parameters)FIEGARCH_L_WW_estimate(method1,wavelet1,level1,parameters,Ydata,ErrorType,zeta);
        parametershat=fmincon(f,parameters0,[],[],[],[],LB,UB,[],options);
        
        % Store estimates in vectors
        d_vec_WW1(i)=parametershat(1);
        omega_vec_WW1(i)=parametershat(2);
        a2_vec_WW1(i)=parametershat(3);       
        b1_vec_WW1(i)=parametershat(4);              
        theta_vec_WW1(i)=parametershat(5);   
        delta_vec_WW1(i)=parametershat(6);  
        
        % Generate vectors for fixed parameters
        if strcmp(ErrorType,'stdged')
            nu_vec_WW1(i)=nu;
        elseif strcmp(ErrorType,'stdstudent')
            nu_vec_WW1(i)=nu;
        elseif strcmp(ErrorType,'skewstudent')
            nu_vec_WW1(i)=nu;
            skew_vec_WW1(i)=skew;
        end;
       
    
    
        %Measure time
        if i<=Measure
         diary Diary_file %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
         diary on;
         disp(datestr(clock,0));
         disp('End Wavelet Whittle1, Start Wavelet Whittle2');
         toc
         tic
         diary off;
        end
     
 %{
    % Wavelet Whittle 2
        % Minimize the Wavelet Whittle objective function
        f=@(parameters)FIEGARCH_L_WW_estimate(method2,wavelet2,level2,parameters,Ydata,ErrorType,zeta);
        parametershat=fmincon(f,parameters0,[],[],[],[],LB,UB,[],options);
        
        % Store estimates in vectors
        d_vec_WW2(i)=parametershat(1);
        omega_vec_WW2(i)=parametershat(2);
        a2_vec_WW2(i)=parametershat(3);       
        b1_vec_WW2(i)=parametershat(4);              
        theta_vec_WW2(i)=parametershat(5);   
        delta_vec_WW2(i)=parametershat(6);  
        
        % Generate vectors for fixed parameters
        if strcmp(ErrorType,'stdged')
            nu_vec_WW2(i)=nu;
        elseif strcmp(ErrorType,'stdstudent')
            nu_vec_WW2(i)=nu;
        elseif strcmp(ErrorType,'skewstudent')
            nu_vec_WW2(i)=nu;
            skew_vec_WW2(i)=skew;
        end;
        
  
        %Measure time
        if i<=Measure
         diary Diary_file %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
         diary on;
         disp(datestr(clock,0));
         disp('End Wavelet Whittle2, Start Wavelet Whittle3');
         toc
         tic
         diary off;
        end
         %}
    
    %{
    % Wavelet Whittle 3
         % Minimize the Wavelet Whittle objective function
        f=@(parameters)FIEGARCH_L_WW_estimate(method3,wavelet3,level3,parameters,Ydata,ErrorType,zeta);
        parametershat=fmincon(f,parameters0,[],[],[],[],LB,UB,[],options);
        
        % Store estimates in vectors
        d_vec_WW3(i)=parametershat(1);
        omega_vec_WW3(i)=parametershat(2);
        a2_vec_WW3(i)=parametershat(3);       
        b1_vec_WW3(i)=parametershat(4);              
        theta_vec_WW3(i)=parametershat(5);   
        delta_vec_WW3(i)=parametershat(6);  
        
        % Generate vectors for fixed parameters
        if strcmp(ErrorType,'stdged')
            nu_vec_WW3(i)=nu;
        elseif strcmp(ErrorType,'stdstudent')
            nu_vec_WW3(i)=nu;
        elseif strcmp(ErrorType,'skewstudent')
            nu_vec_WW3(i)=nu;
            skew_vec_WW3(i)=skew;
        end;
    
        %Measure time
        if i<=Measure
         diary Diary_file %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
         diary on;
         disp(datestr(clock,0));
         disp('End Wavelet Whittle3');
         toc
         tic
         disp('End of the round');
         diary off;
        end
          %}     
    
    % Go to the next simulation (total number of simulations = M)
         
    end    
    matlabpool close
    
elseif zetaflag==1
    disp('You are estimating zeta ')  
 %Extension not used in the paper
  %for i=1:M ;
    
    matlabpool open % lock or unlock both open here and close at the bottom  
    parfor i=1:M
        
        disp(['i=', int2str(i)]);
        
        % Simulate data
        % A) No Jumps
             %[Ydata,ehtSim,ZtSim]=FIEGARCH_simulate(N,d,omega,aL,bL,theta,delta,ErrorType,zeta); % Nx1
        
        % B) Jumps
        % Adjust jump settings in the function
             [Ydata,ehtSim,ZtSim,Jumps]=FIEGARCH_Jump_simulate(N,d,omega,aL,bL,theta,delta,ErrorType,zeta); % Nx1
        
        tic
%{
     % MLE
        %Measure time
        if i<=Measure
         diary Diary_file %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
         diary on;
         disp(datestr(clock,0));
         disp('Start MLE');
         toc
         tic
         diary off;
        end
        
       
        % Minimize negative log-likelihood
        f_MLE=@(parameters)FIEGARCH_L_ML_estimate(parameters,Ydata,ErrorType,zeta);
        parametershat_MLE=fmincon(f_MLE,parameters0,[],[],[],[],LBMLE,UBMLE,[],options);
        
        % Store estimates in vectors
        d_vec_MLE(i)=parametershat_MLE(1);
        omega_vec_MLE(i)=parametershat_MLE(2);
        a2_vec_MLE(i)=parametershat_MLE(3);       
        b1_vec_MLE(i)=parametershat_MLE(4);              
        theta_vec_MLE(i)=parametershat_MLE(5);   
        delta_vec_MLE(i)=parametershat_MLE(6);  
       
        % Generate vectors for fixed parameters
        if strcmp(ErrorType,'stdged')
            nu_vec_MLE(i)=parametershat(7);
        elseif strcmp(ErrorType,'stdstudent')
            nu_vec_MLE(i)=parametershat(7);
        elseif strcmp(ErrorType,'skewstudent')
            nu_vec_MLE(i)=parametershat(7);
            skew_vec_MLE(i)=parametershat(8);
        end;
        
        
        %Measure time
        if i<=Measure
         diary Diary_file %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
         diary on;
         disp('End MLE, Start Whittle');
         disp(datestr(clock,0));
         toc
         tic
         diary off;
        end
         %} 
       
        

    % Whittle
        % Minimize the Whittle objective function
        f_W=@(parameters)FIEGARCH_L_W_estimate(parameters,Ydata,ErrorType,zeta);
        parametershat_W=fmincon(f_W,parameters0,[],[],[],[],LB,UB,[],options);
        
        % Store estimates in vectors
        d_vec_W(i)=parametershat_W(1);
        omega_vec_W(i)=parametershat_W(2);
        a2_vec_W(i)=parametershat_W(3);       
        b1_vec_W(i)=parametershat_W(4);              
        theta_vec_W(i)=parametershat_W(5);   
        delta_vec_W(i)=parametershat_W(6);  
        
        % Generate vectors for fixed parameters
        if strcmp(ErrorType,'stdged')
            nu_vec_W(i)=parametershat_W(7);
        elseif strcmp(ErrorType,'stdstudent')
            nu_vec_W(i)=parametershat_W(7);
        elseif strcmp(ErrorType,'skewstudent')
            nu_vec_W(i)=parametershat_W(7);
            skew_vec_W(i)=parametershat_W(8);
        end;
       
        
        
        %Measure time
        if i<=Measure
         diary Diary_file %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
         diary on;
         disp(datestr(clock,0));
         disp('End Whittle, Start Wavelet Whittle');
         toc
         tic
         diary off;
        end
  
      
    % Wavelet Whittle 1
        % Minimize the Wavelet Whittle objective function
        f=@(parameters)FIEGARCH_L_WW_estimate(method1,wavelet1,level1,parameters,Ydata,ErrorType,zeta);
        parametershat=fmincon(f,parameters0,[],[],[],[],LB,UB,[],options);
        
        % Store estimates in vectors
        d_vec_WW1(i)=parametershat(1);
        omega_vec_WW1(i)=parametershat(2);
        a2_vec_WW1(i)=parametershat(3);       
        b1_vec_WW1(i)=parametershat(4);              
        theta_vec_WW1(i)=parametershat(5);   
        delta_vec_WW1(i)=parametershat(6);  
        
        % Generate vectors for fixed parameters
        if strcmp(ErrorType,'stdged')
            nu_vec_WW1(i)=parametershat(7);
        elseif strcmp(ErrorType,'stdstudent')
            nu_vec_WW1(i)=parametershat(7);
        elseif strcmp(ErrorType,'skewstudent')
            nu_vec_WW1(i)=parametershat(7);
            skew_vec_WW1(i)=parametershat(8);
        end;
       
    
    
        %Measure time
        if i<=Measure
         diary Diary_file %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
         diary on;
         disp(datestr(clock,0));
         disp('End Wavelet Whittle1, Start Wavelet Whittle2');
         toc
         tic
         diary off;
        end
      
  %{
    % Wavelet Whittle 2
        % Minimize the Wavelet Whittle objective function
        f=@(parameters)FIEGARCH_L_WW_estimate(method2,wavelet2,level2,parameters,Ydata,ErrorType,zeta);
        parametershat=fmincon(f,parameters0,[],[],[],[],LB,UB,[],options);
        
        % Store estimates in vectors
        d_vec_WW2(i)=parametershat(1);
        omega_vec_WW2(i)=parametershat(2);
        a2_vec_WW2(i)=parametershat(3);       
        b1_vec_WW2(i)=parametershat(4);              
        theta_vec_WW2(i)=parametershat(5);   
        delta_vec_WW2(i)=parametershat(6);  
        
        % Generate vectors for fixed parameters
        if strcmp(ErrorType,'stdged')
            nu_vec_WW2(i)=parametershat(7);
        elseif strcmp(ErrorType,'stdstudent')
            nu_vec_WW2(i)=parametershat(7);
        elseif strcmp(ErrorType,'skewstudent')
            nu_vec_WW2(i)=parametershat(7);
            skew_vec_WW2(i)=parametershat(8);
        end;
        
  
        %Measure time
        if i<=Measure
         diary Diary_file %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
         diary on;
         disp(datestr(clock,0));
         disp('End Wavelet Whittle2, Start Wavelet Whittle3');
         toc
         tic
         diary off;
        end
         %}
    
    %{
    % Wavelet Whittle 3
         % Minimize the Wavelet Whittle objective function
        f=@(parameters)FIEGARCH_L_WW_estimate(method3,wavelet3,level3,parameters,Ydata,ErrorType,zeta);
        parametershat=fmincon(f,parameters0,[],[],[],[],LB,UB,[],options);
        
        % Store estimates in vectors
        d_vec_WW3(i)=parametershat(1);
        omega_vec_WW3(i)=parametershat(2);
        a2_vec_WW3(i)=parametershat(3);       
        b1_vec_WW3(i)=parametershat(4);              
        theta_vec_WW3(i)=parametershat(5);   
        delta_vec_WW3(i)=parametershat(6);  
        
        % Generate vectors for fixed parameters
        if strcmp(ErrorType,'stdged')
            nu_vec_WW3(i)=parametershat(7);
        elseif strcmp(ErrorType,'stdstudent')
            nu_vec_WW3(i)=parametershat(7);
        elseif strcmp(ErrorType,'skewstudent')
            nu_vec_WW3(i)=parametershat(7);
            skew_vec_WW3(i)=parametershat(8);
        end;
    
        %Measure time
        if i<=Measure
         diary Diary_file %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
         diary on;
         disp(datestr(clock,0));
         disp('End Wavelet Whittle3');
         toc
         tic
         disp('End of the round');
         diary off;
        end
          %}     
    
    % Go to the next simulation (total number of simulations = M)
         
    end    
    matlabpool close

 
end

%MLE
dhat_MLE=mean(d_vec_MLE);
bias1_MLE=dhat_MLE-d;                    
RMSE1_MLE=sqrt(mean((d_vec_MLE-d).^2));  

omegahat_MLE=mean(omega_vec_MLE);
bias2_MLE=omegahat_MLE-omega;                    
RMSE2_MLE=sqrt(mean((omega_vec_MLE-omega).^2));  

a2hat_MLE=mean(a2_vec_MLE);                
bias3_MLE=a2hat_MLE-a2;                    
RMSE3_MLE=sqrt(mean((a2_vec_MLE-a2).^2));  

b1hat_MLE=mean(b1_vec_MLE);                
bias4_MLE=b1hat_MLE-b1;                   
RMSE4_MLE=sqrt(mean((b1_vec_MLE-b1).^2));  

thetahat_MLE=mean(theta_vec_MLE);               
bias5_MLE=thetahat_MLE-theta;                  
RMSE5_MLE=sqrt(mean((theta_vec_MLE-theta).^2)); 

deltahat_MLE=mean(delta_vec_MLE);              
bias6_MLE=deltahat_MLE-delta;                  
RMSE6_MLE=sqrt(mean((delta_vec_MLE-delta).^2)); 

% Whittle
dhat_W=mean(d_vec_W);
bias1_W=dhat_W-d;                    
RMSE1_W=sqrt(mean((d_vec_W-d).^2));  

omegahat_W=mean(omega_vec_W);
bias2_W=omegahat_W-omega;                    
RMSE2_W=sqrt(mean((omega_vec_W-omega).^2));  

a2hat_W=mean(a2_vec_W);                
bias3_W=a2hat_W-a2;                    
RMSE3_W=sqrt(mean((a2_vec_W-a2).^2));  

b1hat_W=mean(b1_vec_W);                
bias4_W=b1hat_W-b1;                   
RMSE4_W=sqrt(mean((b1_vec_W-b1).^2));  

thetahat_W=mean(theta_vec_W);               
bias5_W=thetahat_W-theta;                  
RMSE5_W=sqrt(mean((theta_vec_W-theta).^2)); 

deltahat_W=mean(delta_vec_W);              
bias6_W=deltahat_W-delta;                  
RMSE6_W=sqrt(mean((delta_vec_W-delta).^2));

%Wavelet Whittle1
dhat_WW1=mean(d_vec_WW1);
bias1_WW1=dhat_WW1-d;                    
RMSE1_WW1=sqrt(mean((d_vec_WW1-d).^2));  

omegahat_WW1=mean(omega_vec_WW1);
bias2_WW1=omegahat_WW1-omega;                    
RMSE2_WW1=sqrt(mean((omega_vec_WW1-omega).^2));  

a2hat_WW1=mean(a2_vec_WW1);                
bias3_WW1=a2hat_WW1-a2;                    
RMSE3_WW1=sqrt(mean((a2_vec_WW1-a2).^2));  

b1hat_WW1=mean(b1_vec_WW1);                
bias4_WW1=b1hat_WW1-b1;                   
RMSE4_WW1=sqrt(mean((b1_vec_WW1-b1).^2));  

thetahat_WW1=mean(theta_vec_WW1);               
bias5_WW1=thetahat_WW1-theta;                  
RMSE5_WW1=sqrt(mean((theta_vec_WW1-theta).^2)); 

deltahat_WW1=mean(delta_vec_WW1);              
bias6_WW1=deltahat_WW1-delta;                  
RMSE6_WW1=sqrt(mean((delta_vec_WW1-delta).^2)); 

%Wavelet Whittle2
dhat_WW2=mean(d_vec_WW2);
bias1_WW2=dhat_WW2-d;                    
RMSE1_WW2=sqrt(mean((d_vec_WW2-d).^2));  

omegahat_WW2=mean(omega_vec_WW2);
bias2_WW2=omegahat_WW2-omega;                    
RMSE2_WW2=sqrt(mean((omega_vec_WW2-omega).^2));  

a2hat_WW2=mean(a2_vec_WW2);                
bias3_WW2=a2hat_WW2-a2;                    
RMSE3_WW2=sqrt(mean((a2_vec_WW2-a2).^2));  

b1hat_WW2=mean(b1_vec_WW2);                
bias4_WW2=b1hat_WW2-b1;                   
RMSE4_WW2=sqrt(mean((b1_vec_WW2-b1).^2));  

thetahat_WW2=mean(theta_vec_WW2);               
bias5_WW2=thetahat_WW2-theta;                  
RMSE5_WW2=sqrt(mean((theta_vec_WW2-theta).^2)); 

deltahat_WW2=mean(delta_vec_WW2);              
bias6_WW2=deltahat_WW2-delta;                  
RMSE6_WW2=sqrt(mean((delta_vec_WW2-delta).^2)); 

%Wavelet Whittle3
dhat_WW3=mean(d_vec_WW3);
bias1_WW3=dhat_WW3-d;                    
RMSE1_WW3=sqrt(mean((d_vec_WW3-d).^2));  

omegahat_WW3=mean(omega_vec_WW3);
bias2_WW3=omegahat_WW3-omega;                    
RMSE2_WW3=sqrt(mean((omega_vec_WW3-omega).^2));  

a2hat_WW3=mean(a2_vec_WW3);                
bias3_WW3=a2hat_WW3-a2;                    
RMSE3_WW3=sqrt(mean((a2_vec_WW3-a2).^2));  

b1hat_WW3=mean(b1_vec_WW3);                
bias4_WW3=b1hat_WW3-b1;                   
RMSE4_WW3=sqrt(mean((b1_vec_WW3-b1).^2));  

thetahat_WW3=mean(theta_vec_WW3);               
bias5_WW3=thetahat_WW3-theta;                  
RMSE5_WW3=sqrt(mean((theta_vec_WW3-theta).^2)); 

deltahat_WW3=mean(delta_vec_WW3);              
bias6_WW3=deltahat_WW3-delta;                  
RMSE6_WW3=sqrt(mean((delta_vec_WW3-delta).^2)); 

 %Summary
if strcmp(ErrorType,'stdnormal')
    Row_names={'\hat{d}' '\hat{\omega}' '\hat{\alpha_{2}}' '\hat{\phi_{1}}' '\hat{\theta}' '\hat{\gamma}'};
    %MLE
    parametershat_MLE=[dhat_MLE omegahat_MLE a2hat_MLE b1hat_MLE thetahat_MLE deltahat_MLE ];
    bias_MLE=[bias1_MLE bias2_MLE bias3_MLE bias4_MLE bias5_MLE bias6_MLE];
    RMSE_MLE=[RMSE1_MLE RMSE2_MLE RMSE3_MLE RMSE4_MLE RMSE5_MLE RMSE6_MLE ];
    paramvec_MLE=[d_vec_MLE omega_vec_MLE a2_vec_MLE b1_vec_MLE theta_vec_MLE delta_vec_MLE ];
    nu=[];
    
    
    %Whittle
    parametershat_W=[dhat_W omegahat_W a2hat_W b1hat_W thetahat_W deltahat_W ];
    bias_W=[bias1_W bias2_W bias3_W bias4_W bias5_W bias6_W];
    RMSE_W=[RMSE1_W RMSE2_W RMSE3_W RMSE4_W RMSE5_W RMSE6_W ];
    paramvec_W=[d_vec_W omega_vec_W a2_vec_W b1_vec_W theta_vec_W delta_vec_W ];
    nu=[];

    %Wavelet Whittle1
    parametershat_WW1=[dhat_WW1 omegahat_WW1 a2hat_WW1 b1hat_WW1 thetahat_WW1 deltahat_WW1 ];
    bias_WW1=[bias1_WW1 bias2_WW1 bias3_WW1 bias4_WW1 bias5_WW1 bias6_WW1];
    RMSE_WW1=[RMSE1_WW1 RMSE2_WW1 RMSE3_WW1 RMSE4_WW1 RMSE5_WW1 RMSE6_WW1 ];
    paramvec_WW1=[d_vec_WW1 omega_vec_WW1 a2_vec_WW1 b1_vec_WW1 theta_vec_WW1 delta_vec_WW1 ];
    nu=[];
    
    %Wavelet Whittle2
    parametershat_WW2=[dhat_WW2 omegahat_WW2 a2hat_WW2 b1hat_WW2 thetahat_WW2 deltahat_WW2];
    bias_WW2=[bias1_WW2 bias2_WW2 bias3_WW2 bias4_WW2 bias5_WW2 bias6_WW2];
    RMSE_WW2=[RMSE1_WW2 RMSE2_WW2 RMSE3_WW2 RMSE4_WW2 RMSE5_WW2 RMSE6_WW2];
    paramvec_WW2=[d_vec_WW2 omega_vec_WW2 a2_vec_WW2 b1_vec_WW2 theta_vec_WW2 delta_vec_WW2];

    %Wavelet Whittle3
    parametershat_WW3=[dhat_WW3 omegahat_WW3 a2hat_WW3 b1hat_WW3 thetahat_WW3 deltahat_WW3 ];
    bias_WW3=[bias1_WW3 bias2_WW3 bias3_WW3 bias4_WW3 bias5_WW3 bias6_WW3];
    RMSE_WW3=[RMSE1_WW3 RMSE2_WW3 RMSE3_WW3 RMSE4_WW3 RMSE5_WW3 RMSE6_WW3 ];
    paramvec_WW3=[d_vec_WW3 omega_vec_WW3 a2_vec_WW3 b1_vec_WW3 theta_vec_WW3 delta_vec_WW3 ];
    nu=[];
    
    
elseif strcmp(ErrorType,'stdged')
    Row_names={'\hat{d}}' '\hat{omega}' '\hat{a2}' '\hat{b1}' '\hat{theta}' '\hat{delta}' '\hat{nu}'};
    %MLE
    nuhat_MLE=mean(nu_vec_MLE) ;             
    bias7_MLE=nuhat_MLE-nu;                    
    RMSE7_MLE=sqrt(mean((nu_vec_MLE-nu).^2));
    
    parametershat_MLE=[dhat_MLE omegahat_MLE a2hat_MLE b1hat_MLE thetahat_MLE deltahat_MLE nuhat_MLE];
    bias_MLE=[bias1_MLE bias2_MLE bias3_MLE bias4_MLE bias5_MLE bias6_MLE bias7_MLE];
    RMSE_MLE=[RMSE1_MLE RMSE2_MLE RMSE3_MLE RMSE4_MLE RMSE5_MLE RMSE6_MLE RMSE7_MLE ];
    paramvec_MLE=[d_vec_MLE omega_vec_MLE a2_vec_MLE b1_vec_MLE theta_vec_MLE delta_vec_MLE nu_vec_MLE];
    
    %Whittle
    nuhat_W=mean(nu_vec_W) ;             
    bias7_W=nuhat_W-nu;                    
    RMSE7_W=sqrt(mean((nu_vec_W-nu).^2));
    
    parametershat_W=[dhat_W omegahat_W a2hat_W b1hat_W thetahat_W deltahat_W nuhat_W];
    bias_W=[bias1_W bias2_W bias3_W bias4_W bias5_W bias6_W bias7_W];
    RMSE_W=[RMSE1_W RMSE2_W RMSE3_W RMSE4_W RMSE5_W RMSE6_W RMSE7_W ];
    paramvec_W=[d_vec_W omega_vec_W a2_vec_W b1_vec_W theta_vec_W delta_vec_W nu_vec_W];

    %Wavelet Whittle1
    nuhat_WW1=mean(nu_vec_WW1) ;             
    bias7_WW1=nuhat_WW1-nu;                    
    RMSE7_WW1=sqrt(mean((nu_vec_WW1-nu).^2));
    
    parametershat_WW1=[dhat_WW1 omegahat_WW1 a2hat_WW1 b1hat_WW1 thetahat_WW1 deltahat_WW1 nuhat_WW1];
    bias_WW1=[bias1_WW1 bias2_WW1 bias3_WW1 bias4_WW1 bias5_WW1 bias6_WW1 bias7_WW1];
    RMSE_WW1=[RMSE1_WW1 RMSE2_WW1 RMSE3_WW1 RMSE4_WW1 RMSE5_WW1 RMSE6_WW1 RMSE7_WW1 ];
    paramvec_WW1=[d_vec_WW1 omega_vec_WW1 a2_vec_WW1 b1_vec_WW1 theta_vec_WW1 delta_vec_WW1 nu_vec_WW1];
    
        %Wavelet Whittle2
    nuhat_WW2=mean(nu_vec_WW2);              
    bias7_WW2=nuhat_WW2-nu;                    
    RMSE7_WW2=sqrt(mean((nu_vec_WW2-nu).^2)); 
    
    parametershat_WW2=[dhat_WW2 omegahat_WW2 a2hat_WW2 b1hat_WW2 thetahat_WW2 deltahat_WW2 nuhat_WW2];
    bias_WW2=[bias1_WW2 bias2_WW2 bias3_WW2 bias4_WW2 bias5_WW2 bias6_WW2 bias7_WW2];
    RMSE_WW2=[RMSE1_WW2 RMSE2_WW2 RMSE3_WW2 RMSE4_WW2 RMSE5_WW2 RMSE6_WW2 RMSE7_WW2 ];
    paramvec_WW2=[d_vec_WW2 omega_vec_WW2 a2_vec_WW2 b1_vec_WW2 theta_vec_WW2 delta_vec_WW2 nu_vec_WW2];
   
    %Wavelet Whittle3
    nuhat_WW3=mean(nu_vec_WW3) ;             
    bias7_WW3=nuhat_WW3-nu;                    
    RMSE7_WW3=sqrt(mean((nu_vec_WW3-nu).^2));
    
    parametershat_WW3=[dhat_WW3 omegahat_WW3 a2hat_WW3 b1hat_WW3 thetahat_WW3 deltahat_WW3 nuhat_WW3];
    bias_WW3=[bias1_WW3 bias2_WW3 bias3_WW3 bias4_WW3 bias5_WW3 bias6_WW3 bias7_WW3];
    RMSE_WW3=[RMSE1_WW3 RMSE2_WW3 RMSE3_WW3 RMSE4_WW3 RMSE5_WW3 RMSE6_WW3 RMSE7_WW3 ];
    paramvec_WW3=[d_vec_WW3 omega_vec_WW3 a2_vec_WW3 b1_vec_WW3 theta_vec_WW3 delta_vec_WW3 nu_vec_WW3];

elseif strcmp(ErrorType,'stdstudent') 
    Row_names={'\hat{d}}' '\hat{omega}' '\hat{a2}' '\hat{b1}' '\hat{theta}' '\hat{delta}' '\hat{nu}'};
    %MLE
    nuhat_MLE=mean(nu_vec_MLE);              
    bias7_MLE=nuhat_MLE-nu;                    
    RMSE7_MLE=sqrt(mean((nu_vec_MLE-nu).^2)); 
    
    parametershat_MLE=[dhat_MLE omegahat_MLE a2hat_MLE b1hat_MLE thetahat_MLE deltahat_MLE nuhat_MLE];
    bias_MLE=[bias1_MLE bias2_MLE bias3_MLE bias4_MLE bias5_MLE bias6_MLE bias7_MLE];
    RMSE_MLE=[RMSE1_MLE RMSE2_MLE RMSE3_MLE RMSE4_MLE RMSE5_MLE RMSE6_MLE RMSE7_MLE ];
    paramvec_MLE=[d_vec_MLE omega_vec_MLE a2_vec_MLE b1_vec_MLE theta_vec_MLE delta_vec_MLE nu_vec_MLE];
    
    %Whittle
    nuhat_W=mean(nu_vec_W);              
    bias7_W=nuhat_W-nu;                    
    RMSE7_W=sqrt(mean((nu_vec_W-nu).^2)); 
    
    parametershat_W=[dhat_W omegahat_W a2hat_W b1hat_W thetahat_W deltahat_W nuhat_W];
    bias_W=[bias1_W bias2_W bias3_W bias4_W bias5_W bias6_W bias7_W];
    RMSE_W=[RMSE1_W RMSE2_W RMSE3_W RMSE4_W RMSE5_W RMSE6_W RMSE7_W ];
    paramvec_W=[d_vec_W omega_vec_W a2_vec_W b1_vec_W theta_vec_W delta_vec_W nu_vec_W];
    
    %Wavelet Whittle1
    nuhat_WW1=mean(nu_vec_WW1);              
    bias7_WW1=nuhat_WW1-nu;                    
    RMSE7_WW1=sqrt(mean((nu_vec_WW1-nu).^2)); 
    
    parametershat_WW1=[dhat_WW1 omegahat_WW1 a2hat_WW1 b1hat_WW1 thetahat_WW1 deltahat_WW1 nuhat_WW1];
    bias_WW1=[bias1_WW1 bias2_WW1 bias3_WW1 bias4_WW1 bias5_WW1 bias6_WW1 bias7_WW1];
    RMSE_WW1=[RMSE1_WW1 RMSE2_WW1 RMSE3_WW1 RMSE4_WW1 RMSE5_WW1 RMSE6_WW1 RMSE7_WW1 ];
    paramvec_WW1=[d_vec_WW1 omega_vec_WW1 a2_vec_WW1 b1_vec_WW1 theta_vec_WW1 delta_vec_WW1 nu_vec_WW1];
    
    %Wavelet Whittle2
     nuhat_WW2=mean(nu_vec_WW2);               
    bias7_WW2=nuhat_WW2-nu;                    
    RMSE7_WW2=sqrt(mean((nu_vec_WW2-nu).^2)); 
    
    skewhat_WW=mean(skew_vec_WW2);               
    bias8_WW=skewhat_WW-skew;                    
    RMSE8_WW=sqrt(mean((skew_vec_WW2-skew).^2));
    
    parametershat_WW2=[dhat_WW2 omegahat_WW2 a2hat_WW2 b1hat_WW2 thetahat_WW2 deltahat_WW2 nuhat_WW2 skewhat_WW]
    bias_WW2=[bias1_WW2 bias2_WW2 bias3_WW2 bias4_WW2 bias5_WW2 bias6_WW2 bias7_WW2 ];
    RMSE_WW2=[RMSE1_WW2 RMSE2_WW2 RMSE3_WW2 RMSE4_WW2 RMSE5_WW2 RMSE6_WW2 RMSE7_WW2];
    paramvec_WW2=[d_vec_WW2 omega_vec_WW2 a2_vec_WW2 b1_vec_WW2 theta_vec_WW2 delta_vec_WW2 nu_vec_WW2 skew_vec_WW2];

    %Wavelet Whittle3
    nuhat_WW3=mean(nu_vec_WW3);              
    bias7_WW3=nuhat_WW3-nu;                    
    RMSE7_WW3=sqrt(mean((nu_vec_WW3-nu).^2)); 
    
    parametershat_WW3=[dhat_WW3 omegahat_WW3 a2hat_WW3 b1hat_WW3 thetahat_WW3 deltahat_WW3 nuhat_WW3];
    bias_WW3=[bias1_WW3 bias2_WW3 bias3_WW3 bias4_WW3 bias5_WW3 bias6_WW3 bias7_WW3];
    RMSE_WW3=[RMSE1_WW3 RMSE2_WW3 RMSE3_WW3 RMSE4_WW3 RMSE5_WW3 RMSE6_WW3 RMSE7_WW3 ];
    paramvec_WW3=[d_vec_WW3 omega_vec_WW3 a2_vec_WW3 b1_vec_WW3 theta_vec_WW3 delta_vec_WW3 nu_vec_WW3];
   
    
elseif strcmp(ErrorType,'skewstudent')
    Row_names={'\hat{d}}' '\hat{omega}' '\hat{a2}' '\hat{b1}' '\hat{theta}' '\hat{delta}' '\hat{nu}' '\hat{skew}'};
     %MLE
    nuhat_MLE=mean(nu_vec_MLE);               
    bias7_MLE=nuhat_MLE-nu;                    
    RMSE7_MLE=sqrt(mean((nu_vec_MLE-nu).^2)); 
    
    skewhat_MLE=mean(skew_vec_MLE);               
    bias8_MLE=skewhat_MLE-skew;                    
    RMSE8_MLE=sqrt(mean((skew_vec_MLE-skew).^2));
    
    parametershat_MLE=[dhat_MLE omegahat_MLE a2hat_MLE b1hat_MLE thetahat_MLE deltahat_MLE nuhat_MLE skewhat_MLE]
    bias_MLE=[bias1_MLE bias2_MLE bias3_MLE bias4_MLE bias5_MLE bias6_MLE bias7_MLE bias8_MLE];
    RMSE_MLE=[RMSE1_MLE RMSE2_MLE RMSE3_MLE RMSE4_MLE RMSE5_MLE RMSE6_MLE RMSE7_MLE RMSE8_MLE];
    paramvec_MLE=[d_vec_MLE omega_vec_MLE a2_vec_MLE b1_vec_MLE theta_vec_MLE delta_vec_MLE nu_vec_MLE skew_vec_MLE];
    
    %Whittle
    nuhat_W=mean(nu_vec_W);               
    bias7_W=nuhat_W-nu;                    
    RMSE7_W=sqrt(mean((nu_vec_W-nu).^2)); 
    
    skewhat_W=mean(skew_vec_W);               
    bias8_W=skewhat_W-skew;                    
    RMSE8_W=sqrt(mean((skew_vec_W-skew).^2));
    
    parametershat_W=[dhat_W omegahat_W a2hat_W b1hat_W thetahat_W deltahat_W nuhat_W skewhat_W]
    bias_W=[bias1_W bias2_W bias3_W bias4_W bias5_W bias6_W bias7_W bias8_W];
    RMSE_W=[RMSE1_W RMSE2_W RMSE3_W RMSE4_W RMSE5_W RMSE6_W RMSE7_W RMSE8_W];
    paramvec_W=[d_vec_W omega_vec_W a2_vec_W b1_vec_W theta_vec_W delta_vec_W nu_vec_W skew_vec_W];
    
    %Wavelet Whittle
     nuhat_WW1=mean(nu_vec_WW1);               
    bias7_WW1=nuhat_WW1-nu;                    
    RMSE7_WW1=sqrt(mean((nu_vec_WW1-nu).^2)); 
    
    skewhat_WW1=mean(skew_vec_WW1);               
    bias8_WW1=skewhat_WW1-skew;                    
    RMSE8_WW1=sqrt(mean((skew_vec_WW1-skew).^2));
    
    parametershat_WW1=[dhat_WW1 omegahat_WW1 a2hat_WW1 b1hat_WW1 thetahat_WW1 deltahat_WW1 nuhat_WW1 skewhat_WW1]
    bias_WW1=[bias1_WW1 bias2_WW1 bias3_WW1 bias4_WW1 bias5_WW1 bias6_WW1 bias7_WW1 bias8_WW1];
    RMSE_WW1=[RMSE1_WW1 RMSE2_WW1 RMSE3_WW1 RMSE4_WW1 RMSE5_WW1 RMSE6_WW1 RMSE7_WW1 RMSE8_WW1];
    paramvec_WW1=[d_vec_WW1 omega_vec_WW1 a2_vec_WW1 b1_vec_WW1 theta_vec_WW1 delta_vec_WW1 nu_vec_WW1 skew_vec_WW1];

    %Wavelet Whittle2
     nuhat_WW2=mean(nu_vec_WW2);               
    bias7_WW2=nuhat_WW2-nu;                    
    RMSE7_WW2=sqrt(mean((nu_vec_WW2-nu).^2)); 
    
    skewhat_WW=mean(skew_vec_WW2);               
    bias8_WW2=skewhat_WW-skew;                    
    RMSE8_WW2=sqrt(mean((skew_vec_WW2-skew).^2));
    
    parametershat_WW2=[dhat_WW2 omegahat_WW2 a2hat_WW2 b1hat_WW2 thetahat_WW2 deltahat_WW2 nuhat_WW2 skewhat_WW]
    bias_WW2=[bias1_WW2 bias2_WW2 bias3_WW2 bias4_WW2 bias5_WW2 bias6_WW2 bias7_WW2 bias8_WW2];
    RMSE_WW2=[RMSE1_WW2 RMSE2_WW2 RMSE3_WW2 RMSE4_WW2 RMSE5_WW2 RMSE6_WW2 RMSE7_WW2 RMSE8_WW2];
    paramvec_WW2=[d_vec_WW2 omega_vec_WW2 a2_vec_WW2 b1_vec_WW2 theta_vec_WW2 delta_vec_WW2 nu_vec_WW2 skew_vec_WW2];

    %Wavelet Whittle3
     nuhat_WW3=mean(nu_vec_WW3);               
    bias7_WW3=nuhat_WW3-nu;                    
    RMSE7_WW3=sqrt(mean((nu_vec_WW3-nu).^2)); 
    
    skewhat_WW3=mean(skew_vec_WW3);               
    bias8_WW3=skewhat_WW3-skew;                    
    RMSE8_WW3=sqrt(mean((skew_vec_WW3-skew).^2));
    
    parametershat_WW3=[dhat_WW3 omegahat_WW3 a2hat_WW3 b1hat_WW3 thetahat_WW3 deltahat_WW3 nuhat_WW3 skewhat_WW3]
    bias_WW3=[bias1_WW3 bias2_WW3 bias3_WW3 bias4_WW3 bias5_WW3 bias6_WW3 bias7_WW3 bias8_WW3];
    RMSE_WW3=[RMSE1_WW3 RMSE2_WW3 RMSE3_WW3 RMSE4_WW3 RMSE5_WW3 RMSE6_WW3 RMSE7_WW3 RMSE8_WW3];
    paramvec_WW3=[d_vec_WW3 omega_vec_WW3 a2_vec_WW3 b1_vec_WW3 theta_vec_WW3 delta_vec_WW3 nu_vec_WW3 skew_vec_WW3];

    
end
        
%%%%%%%%%%%%%%
%    DIARY   %
%%%%%%%%%%%%%%

diary Diary_file %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
diary on;
 disp(datestr(clock,0));
 level1
 level2
 level3
 disp('PARAMETERS, ESTIMATED PARAMETERS, BIAS, RMSE'); 
 disp('******************************************************************');
 disp(['N=',num2str(N),'; M=',num2str(M),'d=', num2str(d),'; omega=',num2str(omega),'; a2=',num2str(a2),'; b1=',num2str(b1),'; theta=',num2str(theta),'; delta=',num2str(delta)]);
 disp(['method1=',method1,'; method2=',method2,'; method3=',method3,'; wavelet1=',wavelet1,'; wavelet2=',wavelet2,'; wavelet3=',wavelet3,'level1=', num2str(level1),'; level2=',num2str(level2),'; level3=',num2str(level3)]);
 disp('******************************************************************');
 disp('RESTRICTIONS -0.5<d<0.5 -inf<omega<inf -1<a1<1,-1<b1<1,0<theta2<inf,-inf<delta<inf');
 disp('  ');
 diary off;
 
if strcmp(ErrorType,'stdnormal')
   parameters=[d,omega,aL,bL,theta,delta];
 diary Diary_file %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
 diary on;
 disp('parameters = d        omega    a2       b1       theta    delta');
 disp('  ');
 disp(['param_true= ',num2str(parameters,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp('  ');
 disp(['param_MLE = ',num2str(parametershat_MLE,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp(['bias_MLE =  ',num2str(bias_MLE,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp(['RMSE_MLE =  ',num2str(RMSE_MLE,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp('  ');
 disp(['param_W =   ',num2str(parametershat_W,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp(['bias_W =    ',num2str(bias_W,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp(['RMSE_W =    ',num2str(RMSE_W,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp('  ');
 disp(['param_hat_WW1 = ',num2str(parametershat_WW1,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp(['bias_WW1 =      ',num2str(bias_WW1,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp(['RMSE_WW1 =      ',num2str(RMSE_WW1,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp('  ');
  disp(['param_hat_WW2 = ',num2str(parametershat_WW2,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp(['bias_WW2 =      ',num2str(bias_WW2,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp(['RMSE_WW2 =      ',num2str(RMSE_WW2,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp('  ');
  disp(['param_hat_WW3 = ',num2str(parametershat_WW3,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp(['bias_WW3 =      ',num2str(bias_WW3,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp(['RMSE_WW3 =      ',num2str(RMSE_WW3,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp('  ');
 disp([datestr(clock,0),'; Time in seconds elapsed in execution =', num2str(etime(clock,t0))]);
 disp('  ');
 
  disp('******************************************************************');
 disp(['N=',num2str(N),'; M=',num2str(M),'d=', num2str(d),'; omega=',num2str(omega),'; a2=',num2str(a2),'; b1=',num2str(b1),'; theta=',num2str(theta),'; delta=',num2str(delta)]);
 disp('******************************************************************');
 disp(' MLE ');
 disp('d        omega    a2       b1       theta    delta');
 disp(num2str(parametershat_MLE, '%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f'));
 
 disp(' Whittle ');
 disp('d        omega    a2       b1       theta    delta');
 disp(num2str(parametershat_W, '%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f'));
 disp(' Wavelet Whittle1 ');
 disp('d        omega    a2       b1       theta    delta');
 disp(num2str(parametershat_WW1, '%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f'));
  disp(' Wavelet Whittle2 ');
 disp('d        omega    a2       b1       theta    delta');
 disp(num2str(parametershat_WW2, '%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f'));
  disp(' Wavelet Whittle3 ');
 disp('d        omega    a2       b1       theta    delta');
 disp(num2str(parametershat_WW3, '%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f'));
 
 disp('  ');

 diary off;

    
elseif strcmp(ErrorType,'stdged')
 parameters=[d,omega,aL,bL,theta,delta,nu];
 diary Diary_file %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
 diary on;
 disp('parameters = d        omega    a2       b1       theta    delta    (nu)');
 disp(['param_true= ',num2str(parameters,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp('  ');
 disp(['param_MLE = ',num2str(parametershat_MLE,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp(['bias_MLE =  ',num2str(bias_MLE,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp(['RMSE_MLE =  ',num2str(RMSE_MLE,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp('  ');
 disp(['param_W =   ',num2str(parametershat_W,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp(['bias_W =    ',num2str(bias_W,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp(['RMSE_W =    ',num2str(RMSE_W,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp('  ');
 disp(['param_hat_WW1 = ',num2str(parametershat_WW1,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp(['bias_WW1 =      ',num2str(bias_WW1,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp(['RMSE_WW1 =      ',num2str(RMSE_WW1,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp('  ');
 disp(['param_hat_WW2 = ',num2str(parametershat_WW2,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp(['bias_WW2 =      ',num2str(bias_WW2,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp(['RMSE_WW2 =      ',num2str(RMSE_WW2,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp('  ');
  disp(['param_hat_WW3 = ',num2str(parametershat_WW3,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp(['bias_WW3 =      ',num2str(bias_WW3,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp(['RMSE_WW3 =      ',num2str(RMSE_WW3,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp('  ');
 disp([datestr(clock,0),'; Time in seconds elapsed in execution =', num2str(etime(clock,t0))]);
 disp('  ');

 disp('******************************************************************');
 disp(['N=',num2str(N),'; M=',num2str(M),'d=', num2str(d),'; omega=',num2str(omega),'; a2=',num2str(a2),'; b1=',num2str(b1),'; theta=',num2str(theta),'; delta=',num2str(delta),'; nu= ',num2str(nu)]);
 disp('******************************************************************');
 disp(' MLE ');
 disp('d        omega    a2       b1       theta    delta    (nu)');
 disp(num2str(parametershat_MLE, '%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f'));
 
 disp(' Whittle ');
 disp('d        omega    a2       b1       theta    delta    (nu)');
 disp(num2str(parametershat_W, '%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f'));
 
 disp('Wavelet Whittle1  ');
 disp('d        omega    a2       b1       theta    delta    (nu)');
 disp(num2str(parametershat_WW1, '%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f'));
 
  disp(' Wavelet Whittle2 ');
 disp('d        omega    a2       b1       theta    delta    (nu)');
 disp(num2str(parametershat_WW2, '%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f'));
 
  disp('Wavelet Whittle3  ');
 disp('d        omega    a2       b1       theta    delta    (nu)');
 disp(num2str(parametershat_WW3, '%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f'));


 disp('  ');
 diary off;

elseif strcmp(ErrorType,'stdstudent')
    parameters=[d,omega,aL,bL,theta,delta,nu];
 diary Diary_file %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
 diary on;
 disp('parameters = d        omega    a2       b1       theta    delta    (nu)');
 disp(['param_true= ',num2str(parameters,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp('  ');
 disp(['param_MLE = ',num2str(parametershat_MLE,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp(['bias_MLE =  ',num2str(bias_MLE,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp(['RMSE_MLE =  ',num2str(RMSE_MLE,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp('  ');
 disp(['param_W =   ',num2str(parametershat_W,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp(['bias_W =    ',num2str(bias_W,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp(['RMSE_W =    ',num2str(RMSE_W,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp('  ');
 disp(['param_hat_WW1 = ',num2str(parametershat_WW1,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp(['bias_WW1 =      ',num2str(bias_WW1,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp(['RMSE_WW1 =      ',num2str(RMSE_WW1,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp('  ');
  disp(['param_hat_WW2 = ',num2str(parametershat_WW2,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp(['bias_WW2 =      ',num2str(bias_WW2,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp(['RMSE_WW2 =      ',num2str(RMSE_WW2,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp('  ');
  disp(['param_hat_WW3 = ',num2str(parametershat_WW3,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp(['bias_WW3 =      ',num2str(bias_WW3,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp(['RMSE_WW3 =      ',num2str(RMSE_WW3,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp('  ');
 disp([datestr(clock,0),'; Time in seconds elapsed in execution =', num2str(etime(clock,t0))]);
 disp('  ');

 disp('******************************************************************');
 disp(['N=',num2str(N),'; M=',num2str(M),'d=', num2str(d),'; omega=',num2str(omega),'; a2=',num2str(a2),'; b1=',num2str(b1),'; theta=',num2str(theta),'; delta=',num2str(delta),'; nu= ',num2str(nu)]);
 disp('******************************************************************');
 disp(' MLE ');
 disp('d        omega    a2       b1       theta    delta    (nu)');
 disp(num2str(parametershat_MLE, '%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f'));
 
 disp(' Whittle ');
 disp('d        omega    a2       b1       theta    delta    (nu)');
 disp(num2str(parametershat_W, '%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f'));
 
 disp(' Wavelet Whittle1 ');
 disp('d        omega    a2       b1       theta    delta    (nu)');
 disp(num2str(parametershat_WW1, '%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f'));
 
  disp(' Wavelet Whittle2 ');
 disp('d        omega    a2       b1       theta    delta    (nu)');
 disp(num2str(parametershat_WW2, '%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f'));
 
  disp(' Wavelet Whittle3 ');
 disp('d        omega    a2       b1       theta    delta    (nu)');
 disp(num2str(parametershat_WW3, '%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f'));

 disp('  ');
 diary off;
    
elseif strcmp(ErrorType,'skewstudent')
    parameters=[d,omega,aL,bL,theta,delta,nu,skew];
 diary Diary_file %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
 diary on;
 disp('parameters = d        omega    a2       b1       theta    delta    (nu)     (skew)');
 disp(['param_true= ',num2str(parameters,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f ')]);
 disp('  ');
 disp(['param_MLE = ',num2str(parametershat_MLE,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp(['bias_MLE =  ',num2str(bias_MLE,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp(['RMSE_MLE =  ',num2str(RMSE_MLE,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 
 disp(['param_W =   ',num2str(parametershat_W,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp(['bias_W =    ',num2str(bias_W,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp(['RMSE_W =    ',num2str(RMSE_W,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 
 disp(['param_hat_WW1 = ',num2str(parametershat_WW1,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp(['bias_WW1 =      ',num2str(bias_WW1,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp(['RMSE_WW1 =      ',num2str(RMSE_WW1,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp('  ');
  disp(['param_hat_WW2 = ',num2str(parametershat_WW2,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp(['bias_WW2 =      ',num2str(bias_WW2,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp(['RMSE_WW2 =      ',num2str(RMSE_WW2,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp('  ');
  disp(['param_hat_WW3 = ',num2str(parametershat_WW3,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp(['bias_WW3 =      ',num2str(bias_WW3,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp(['RMSE_WW3 =      ',num2str(RMSE_WW3,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f')]);
 disp('  ');
 disp([datestr(clock,0),'; Time in seconds elapsed in execution =', num2str(etime(clock,t0))]);
 disp('  ');
 
 disp('******************************************************************');
 disp(['N=',num2str(N),'; M=',num2str(M),'d=', num2str(d),'; omega=',num2str(omega),'; a2=',num2str(a2),'; b1=',num2str(b1),'; theta=',num2str(theta),'; delta=',num2str(delta),'; nu= ',num2str(nu),';lambda=',num2str(skew)]);
 disp('******************************************************************');
 disp(' MLE ');
 disp('d        omega    a2       b1       theta    delta    (nu)     (skew)');
 disp(num2str(parametershat_MLE, '%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f'));
 disp('  ');
 
 disp(' Whittle ');
 disp('d        omega    a2       b1       theta    delta    (nu)     (skew)');
 disp(num2str(parametershat_W, '%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f'));
 
 disp('Wavelet Whittle1 ');
 disp('d        omega    a2       b1       theta    delta    (nu)     (skew)');
 disp(num2str(parametershat_WW1, '%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f'));
 
  disp('Wavelet Whittle2 ');
 disp('d        omega    a2       b1       theta    delta    (nu)     (skew)');
 disp(num2str(parametershat_WW2, '%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f'));

  disp('Wavelet Whittle3 ');
 disp('d        omega    a2       b1       theta    delta    (nu)     (skew)');
 disp(num2str(parametershat_WW3, '%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f'));

 diary off;
 
end

%Output code: generate for TEX

Column_names={'TRUE' 'METHOD' 'MEAN' 'BIAS' 'RMSE'};
Estimators={'WWE1' 'WWE2' 'WWE3' 'FWE' 'MLE'};

%param=[];
hat=zeros(5*length(parameters),1);
rmse=zeros(5*length(parameters),1);
bias=zeros(5*length(parameters),1);

i=0;
for t=1:(length(parametershat_W));
    i=i+1;
    %param(i)=Row_names{t}
    
    hat(i)=parametershat_WW1(t);
    rmse(i)=RMSE_WW1(t);
    bias(i)=bias_WW1(t);
    
    i=i+1;
   % param(i)='';
   
     hat(i)=parametershat_WW2(t);
    rmse(i)=RMSE_WW2(t);
    bias(i)=bias_WW2(t);
    i=i+1;
    %param(i)='';
    
     hat(i)=parametershat_WW3(t);
    rmse(i)=RMSE_WW3(t);
    bias(i)=bias_WW3(t);
    i=i+1;
    %param(i)='';
    
     hat(i)=parametershat_W(t);
    rmse(i)=RMSE_W(t);
    bias(i)=bias_W(t);
    i=i+1;
   % param(i)='';
  
     hat(i)=parametershat_MLE(t);
    rmse(i)=RMSE_MLE(t);
    bias(i)=bias_MLE(t);
end

diary Diary_file %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
diary on;

final=[ hat bias rmse];

fprintf('%20s & ','Observations');

fprintf('%7s & ',' ');

fprintf('%7s & ',' ');
fprintf('%1s & ',' ');
fprintf('%7s & ',num2str(N));

fprintf('%7s & ',' ');
fprintf('%7s & ',' ');
fprintf('%7s \\\\ \n',' ');


fprintf('%20s & ','');

for m=1:size(Column_names,2)-1
    fprintf('%7s & ',cell2mat(Column_names(m)))
    if m==2
        fprintf('%1s & ',' ');
    end
        
end
fprintf('%7s & ',cell2mat(Column_names(size(Column_names,2))))
fprintf('%7s \\\\ \n',' ')

for ROW=1:size(final, 1)
 
     fprintf('%20s & ',Row_names{floor((ROW-1)/5)+1});
     
     fprintf('%7.3f & ',parameters(floor((ROW-1)/5)+1));
     
     fprintf('%7s & ',Estimators{ROW-(5*floor((ROW-1)/5))});
     fprintf('%1s & ',' ');
    
for COL=1:size(final,2)-1

if COL==size(Column_names,2)+1
fprintf('%1s & ',' ');
end

  fprintf('%7.3f & ', final(ROW,COL));
  
end
fprintf('%7.3f & ', final(ROW,end));
fprintf('%7s \\\\ \n', ' ');

end
disp('WWE_Donoho 0.25, 2048, Jumps 0.028, N(0,0.2)')% Add montecarlo parameters to distinguish between results in the diary (e.g. jumps parameters and WWE threshold)
diary off;


