
function FIEGARCH_montecarlo_forecast_correction1(scenario)
display('FIEGARCH_montecarlo_forecast_correction1(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 and jump detection)+dealing with jumps in all three estimators!!!
%               Jump detection parts of the code... and jump settings
%               in FIEGARCH_Jump_Simulate, choose FIEGARCH_simulate without
%               jumps line 564 and lock jump detection in 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
%               Choose number of in-sample (insample) and out of sample (N_next) forecasts 

% in case of using other than stdnormal distrubutions corrected measures
% for parameters have to be added in montecarlo summary-just some copy-paste
% work... now it works for stdnormal only

%%%%%%%%%%%%
%   USAGE  %
%%%%%%%%%%%%
% Perform MonteCarlo Experiments on the finite sample properties of the
% Maximum likelihood, 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.

% Forecasting performance is studied and the results are corrected to
% prevet inf and NaN values in case of rare diveregence of estimates. All
% these cases are counted to get a good summary of methods performance


%%%%%%%%%%%%%%%%
%   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='DWT';
    method3='DWT';
    
    wavelet1='d4';
    wavelet2='d4';
    wavelet3='la8';
    
    level1=9; % J-3
    level2=9;
    level3=9;
    
    % insample
    insample=N;
    % out of sample
    N_next=1; 
    
    
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='DWT';
    method3='DWT';
    
    wavelet1='d4';
    wavelet2='d4';
    wavelet3='la8';
    
    level1=12;
    level2=12;
    level3=12;
 
    % insample
    insample=N;
    % out of sample
    N_next=1; 
    
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='DWT';
    method3='DWT';
    
    wavelet1='d4';
    wavelet2='d4';
    wavelet3='la8';
    
    level1=9;
    level2=9;
    level3=9;
 
        % insample
    insample=N;
    % out of sample
    N_next=1; 
    
elseif scenario==4
    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.45;
    aL=0.5;
    bL=0.5;
    theta=-0.3;
    delta=0.5;
    
    method1='modwt';
    method2='DWT';
    method3='DWT';
    
    wavelet1='d4';
    wavelet2='d4';
    wavelet3='la8';
    
    level1=12;
    level2=12;
    level3=12;
    
     % insample
    insample=N;
    % out of sample
    N_next=1; 
    
elseif scenario==5
    M=5;    %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='modwt';
    method2='DWT';
    method3='DWT';
    
    wavelet1='d4';
    wavelet2='d4';
    wavelet3='LA8';
    
    level1=9;
    level2=9;
    level3='dwtcut4';
    
       % insample
    insample=N;
    % out of sample
    N_next=1; 
 
elseif scenario==6
    M=10;    %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='modwt';
    method2='modwt';
    method3='modwt';
    
    wavelet1='d4';
    wavelet2='d4';
    wavelet3='Haar';
    
    level1=9;
    level2=6;
    level3=9; 
    
        % insample
    insample=N;
    % out of sample
    N_next=1; 
    
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' and 'modwt' (use 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);

est_err_MLE_in=zeros(M,1);
est_mse_MLE_in=zeros(M,1);
est_mad_MLE_in=zeros(M,1);
est_err_MLE_out=zeros(M,1);
est_mse_MLE_out=zeros(M,1);
est_mad_MLE_out=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);

est_err_W_in=zeros(M,1);
est_mse_W_in=zeros(M,1);
est_mad_W_in=zeros(M,1);
est_err_W_out=zeros(M,1);
est_mse_W_out=zeros(M,1);
est_mad_W_out=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);

est_err_WW1_in=zeros(M,1);
est_mse_WW1_in=zeros(M,1);
est_mad_WW1_in=zeros(M,1);
est_err_WW1_out=zeros(M,1);
est_mse_WW1_out=zeros(M,1);
est_mad_WW1_out=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);

est_err_WW2_in=zeros(M,1);
est_mse_WW2_in=zeros(M,1);
est_mad_WW2_in=zeros(M,1);
est_err_WW2_out=zeros(M,1);
est_mse_WW2_out=zeros(M,1);
est_mad_WW2_out=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);

est_err_WW3_in=zeros(M,1);
est_mse_WW3_in=zeros(M,1);
est_mad_WW3_in=zeros(M,1);
est_err_WW3_out=zeros(M,1);
est_mse_WW3_out=zeros(M,1);
est_mad_WW3_out=zeros(M,1);

NextVals=zeros(M,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  RESTRICTIONS ON ESTIMATES LB,UB  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%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, somewhat more
% reasonable intervals chosen
LB=[ -0.46 0  -0.9 -0.9 -0.9 -0.9]; % Lower bound (d,omega,a2,b1,theta,delta)
UB=[ 0.46 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.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.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,'skewstudent')
                 LB=[LB 2.1 -0.9];    % Lower bounds for  (d,omega,a2,b1,theta,delta,nu,skew)
                 UB=[UB 10 0.9];     % 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- quite narrow band since we can get first
% guess using sample mean...I set true +/- 3, FWE and WWE cannot estimate
% the mean
LBMLE=[ LB(1) -10 LB(3:length(LB))]; % Lower bound (d,omega,a2,b1,theta,delta(,nu,skew))
UBMLE=[ UB(1)  -4 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','GradObj','on','Display','off','TolFun',1e-4); %perez
%MLE
%options=optimset('LargeScale','on','Algorithm','trust-region-reflective','Display','off');
options=optimset('LargeScale','off','Algorithm','active-set','Display','off');%'interior-point'
%options=optimset('LargeScale','on','Display','off','Diagnostics','on'); %active-set when
%gradient not provided

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

diary Diary_file %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
diary on;
disp(['scenario= ',int2str(scenario)]);
disp('Montecarlo: start');
t0=clock;
disp(datestr(clock,0));
disp('  ')
diary off;

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


if zetaflag==0
              JumpsSim='No Jumps';
              %JumpsSim='Jumps (0.028, N(0,0.2))';
              %JumpsSim='Jumps (0.007, N(0,3))';
              
              JumpsEst='Jumps Not Corrected';                            
              %JumpsEst='Jumps Corrected';
   % 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 
        % Detection in estimators! Turn on or off.. but
        % never forget to set what you really want - in case of detection
        % and forecasting the detection has to be done here..
            
     [Ydata_plus,HtSim_plus,~]=FIEGARCH_simulate(N+N_next,d,omega,aL,bL,theta,delta,ErrorType,zeta); % Nx1
     
   %{  
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     % JUMP DETECTION : here on Ydata_plus because of forecasting and omg%
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%v
    wavelet='d4'; % Filters Available:wtfilter('list') 
    %be careful about the number of coefficients the jump influences (Haar:2-symetrically, D4:3 - sum of 2 smaller effects=-bigger effect)
    boundary='circular'; 
    opts='';

    [WJt, VJt, att] = modwt(Ydata_plus,wavelet,1,boundary,opts);  % Modwt transform
    [WJt, VJt] = modwt_cir_shift(WJt, VJt, wavelet, 1); % Circular shift to match true indices
    Detail=WJt(1:size(WJt),1); % Get wavelet coefficients at first level
          
    
  % Simple Jump location detection: Hard thresholding on first level coeffs
        % Donoho and Johnstone 1994 threshold
         thrs=(2*log(length(Detail)))^(0.5)*2^(0.5)*median(abs(Detail))/0.6745;                                                     
          
         jump_loc=abs(Detail)>thrs;
         nojump_loc=abs(jump_loc-1);
               
         Yavg=[(Ydata_plus(end)+Ydata_plus(2))/2;(Ydata_plus(1:end-2)+Ydata_plus(3:end))/2;(Ydata_plus(end-1)+Ydata_plus(1))/2]; % jump correction for all locations
         
         Ydata_plus=nojump_loc.*Ydata_plus+jump_loc.*Yavg; % do correction at jump_loc indices         
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
       %} 
    
              Ydata=Ydata_plus(1:N);
              ehtSim_plus=exp(HtSim_plus);
              NextVals(i)=mean(ehtSim_plus(N+1:end)) ; 
           
        %{   
     
  % B) Jumps 
        % Detection in estimators! Turn on or off.. but
        % never forget to set what you really want
        
         % Adjust jump settings in the function
        
             [Ydata_plus,HtSim_plus,~,~]=FIEGARCH_Jump_simulate(N+N_next,d,omega,aL,bL,theta,delta,ErrorType,zeta); % Nx      
        
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %      JUMP DETECTION : here on Ydata_plus because of forecasting and omg                                     %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    wavelet='d4'; % Filters Available:wtfilter('list') 
    %be careful about the number of coefficients the jump influences (Haar:2-symetrically, D4:3 - sum of 2 smaller effects=-bigger effect)
    boundary='circular'; 
    opts='';

    [WJt, VJt, att] = modwt(Ydata_plus,wavelet,1,boundary,opts);  % Modwt transform
    [WJt, VJt] = modwt_cir_shift(WJt, VJt, wavelet, 1); % Circular shift to match true indices
    Detail=WJt(1:size(WJt),1); % Get wavelet coefficients at first level
           
      % Simple Jump location detection: Hard thresholding on first level coeffs
        % Donoho and Johnstone 1994 threshold
         thrs=(2*log(length(Detail)))^(0.5)*2^(0.5)*median(abs(Detail))/0.6745;                                                     
          
         jump_loc=abs(Detail)>thrs;
         nojump_loc=abs(jump_loc-1);
               
         %Yavg=[(Ydata_plus(end)+Ydata_plus(2))/2;(Ydata_plus(1:end-2)+Ydata_plus(3:end))/2;(Ydata_plus(end-1)+Ydata_plus(1))/2]; % jump correction for all locations
         
         YdataNJ=Ydata_plus.*nojump_loc;
         Yavg=[(YdataNJ(end)+YdataNJ(2))/2;(YdataNJ(1:end-2)+YdataNJ(3:end))/2;(YdataNJ(end-1)+YdataNJ(1))/2]; 
         
         Ydata_plus=nojump_loc.*Ydata_plus+jump_loc.*Yavg; % do correction at jump_loc indices         
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
              
             
              Ydata=Ydata_plus(1:N);
              ehtSim_plus=exp(HtSim_plus);
              NextVals(i)=mean(ehtSim_plus(N+1:end)) ; 
          %}    
              
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
             
              
  % omega estimate for FWE and WWE separately from other estimates 
   omg=mean(log(Ydata.^2))+1.2686;  % ln(yt^2)=ln(zt^2)+omega+inf mean zero MA
        
        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_MLEe=fmincon(f_MLE,parameters0,[],[],[],[],LBMLE,UBMLE,[],options);
        
        
        % Store estimates in vectors
        d_vec_MLE(i)=parametershat_MLEe(1);
        omega_vec_MLE(i)=parametershat_MLEe(2);
        a2_vec_MLE(i)=parametershat_MLEe(3);       
        b1_vec_MLE(i)=parametershat_MLEe(4) ;             
        theta_vec_MLE(i)=parametershat_MLEe(5);   
        delta_vec_MLE(i)=parametershat_MLEe(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;
        
        % Estimate the eht series (conditional variance (not log-variance)
        %of the process) and calculate mean error for given simulation, mean absolute deviation and
        %mean squared error. This will be summed up later on and the overall RMSE will
        %be calculated (choose forecast interval carefully :) )
        
        [~,eht_f,~]=FIEGARCH_core(Ydata_plus,d_vec_MLE(i),omega_vec_MLE(i),a2_vec_MLE(i),b1_vec_MLE(i),theta_vec_MLE(i),delta_vec_MLE(i), ErrorType,zeta);
        
        est_err_MLE=eht_f-ehtSim_plus;  % number, +/-inf, NaN     
        est_err_MLE_in_Q=est_err_MLE(end-insample-N_next+1:end-N_next);
        
        detectnansMLE=sum(isnan(est_err_MLE_in_Q))
        if detectnansMLE>0
         est_err_MLE_in_099=est_err_MLE_in_Q %continue with NaNs- will be
                                             %eliminated later on
        else
         est_err_MLE_in_099=est_err_MLE_in_Q(abs(est_err_MLE_in_Q)<quantile(abs(est_err_MLE_in_Q),0.99));
        end % if not more than x% of errors extreme, may not be considered
            %poor fit: take x% of smallest errors only
        
        est_err_MLE_in(i)=mean(est_err_MLE_in_099);
        est_mad_MLE_in(i)=mean(abs(est_err_MLE_in_099));
        est_mse_MLE_in(i)=mean((est_err_MLE_in_099).^2);
        
        est_err_MLE_out(i)=mean(est_err_MLE(end-N_next+1:end)); % here NaN can still occur - eliminated later on...
        est_mad_MLE_out(i)=mean(abs(est_err_MLE(end-N_next+1:end)));
        est_mse_MLE_out(i)=mean((est_err_MLE(end-N_next+1:end)).^2);
    
        %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                    % see explanation of the procedure at MLE
        % Minimize the Whittle objective function
        f_W=@(parameters)FIEGARCH_L_W_estimate(parameters,Ydata,ErrorType,zeta);
        parametershat_We=fmincon(f_W,parameters0,[],[],[],[],LB,UB,[],options);
        
        % Store estimates in vectors
        d_vec_W(i)=parametershat_We(1);
        omega_vec_W(i)=omg;
        a2_vec_W(i)=parametershat_We(3);       
        b1_vec_W(i)=parametershat_We(4);              
        theta_vec_W(i)=parametershat_We(5);   
        delta_vec_W(i)=parametershat_We(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;


        [~,eht_f,~]=FIEGARCH_core(Ydata_plus,d_vec_W(i),omega_vec_W(i),a2_vec_W(i),b1_vec_W(i),theta_vec_W(i),delta_vec_W(i), ErrorType,zeta);

        est_err_W=eht_f-ehtSim_plus;       
        est_err_W_in_Q=est_err_W(end-insample-N_next+1:end-N_next);
        
        detectnansW=sum(isnan(est_err_W_in_Q))
        if detectnansW>0
         est_err_W_in_099=est_err_W_in_Q
        else
         est_err_W_in_099=est_err_W_in_Q(abs(est_err_W_in_Q)<quantile(abs(est_err_W_in_Q),0.99));
        end
        
        est_err_W_in(i)=mean(est_err_W_in_099);
        est_mad_W_in(i)=mean(abs(est_err_W_in_099));
        est_mse_W_in(i)=mean((est_err_W_in_099).^2);
        
        est_err_W_out(i)=mean(est_err_W(end-N_next+1:end));
        est_mad_W_out(i)=mean(abs(est_err_W(end-N_next+1:end)));
        est_mse_W_out(i)=mean((est_err_W(end-N_next+1:end)).^2);
        
        
        %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                   % see explanation of the procedure at MLE
        % Minimize the Wavelet Whittle objective function
        f=@(parameters)FIEGARCH_L_WW_estimate(method1,wavelet1,level1,parameters,Ydata,ErrorType,zeta);
        parametershat_WW1e=fmincon(f,parameters0,[],[],[],[],LB,UB,[],options);
        
        % Store estimates in vectors
        d_vec_WW1(i)=parametershat_WW1e(1);
        omega_vec_WW1(i)=omg;
        a2_vec_WW1(i)=parametershat_WW1e(3);       
        b1_vec_WW1(i)=parametershat_WW1e(4);              
        theta_vec_WW1(i)=parametershat_WW1e(5);   
        delta_vec_WW1(i)=parametershat_WW1e(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;
       
        [~,eht_f,~]=FIEGARCH_core(Ydata_plus,d_vec_WW1(i),omega_vec_WW1(i),a2_vec_WW1(i),b1_vec_WW1(i),theta_vec_WW1(i),delta_vec_WW1(i), ErrorType,zeta)
        est_err_WW1=eht_f-ehtSim_plus;
        
        est_err_WW1_in_Q=est_err_WW1(end-insample-N_next+1:end-N_next);
        
        detectnansWW1=sum(isnan(est_err_WW1_in_Q))
        if detectnansWW1>0
         est_err_WW1_in_099=est_err_WW1_in_Q
        else
         est_err_WW1_in_099=est_err_WW1_in_Q(abs(est_err_WW1_in_Q)<quantile(abs(est_err_WW1_in_Q),0.99));
        end
      
        est_err_WW1_in(i)=mean(est_err_WW1_in_099);
        est_mad_WW1_in(i)=mean(abs(est_err_WW1_in_099));
        est_mse_WW1_in(i)=mean((est_err_WW1_in_099).^2);
        
        est_err_WW1_out(i)=mean(est_err_WW1(end-N_next+1:end));
        est_mad_WW1_out(i)=mean(abs(est_err_WW1(end-N_next+1:end)));
        est_mse_WW1_out(i)=mean((est_err_WW1(end-N_next+1:end)).^2);
        
        %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                   % see explanation of the procedure at MLE
        % Minimize the Wavelet Whittle objective function
        f=@(parameters)FIEGARCH_L_WW_estimate(method2,wavelet2,level2,parameters,Ydata,ErrorType,zeta);
        parametershat_WW2e=fmincon(f,parameters0,[],[],[],[],LB,UB,[],options);
        
        % Store estimates in vectors
        d_vec_WW2(i)=parametershat_WW2e(1);
        omega_vec_WW2(i)=omg;
        a2_vec_WW2(i)=parametershat_WW2e(3);       
        b1_vec_WW2(i)=parametershat_WW2e(4);              
        theta_vec_WW2(i)=parametershat_WW2e(5);   
        delta_vec_WW2(i)=parametershat_WW2e(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;

        [~,eht_f,~]=FIEGARCH_core(Ydata_plus,d_vec_WW2(i),omega_vec_WW2(i),a2_vec_WW2(i),b1_vec_WW2(i),theta_vec_WW2(i),delta_vec_WW2(i), ErrorType,zeta);

        est_err_WW2=eht_f-ehtSim_plus;
        
        est_err_WW2_in_Q=est_err_WW2(end-insample-N_next+1:end-N_next);
        
         detectnansWW2=sum(isnan(est_err_WW2_in_Q))
        if detectnansWW2>0
         est_err_WW2_in_099=est_err_WW2_in_Q
        else
         est_err_WW2_in_099=est_err_WW2_in_Q(abs(est_err_WW2_in_Q)<quantile(abs(est_err_WW2_in_Q),0.99));
        end
        
        
        est_err_WW2_in(i)=mean(est_err_WW2_in_099);
        est_mad_WW2_in(i)=mean(abs(est_err_WW2_in_099));
        est_mse_WW2_in(i)=mean((est_err_WW2_in_099).^2);
        
        est_err_WW2_out(i)=mean(est_err_WW2(end-N_next+1:end));
        est_mad_WW2_out(i)=mean(abs(est_err_WW2(end-N_next+1:end)));
        est_mse_WW2_out(i)=mean((est_err_WW2(end-N_next+1:end)).^2);
        
        %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                   % see explanation of the procedure at MLE
         % Minimize the Wavelet Whittle objective function
        f=@(parameters)FIEGARCH_L_WW_estimate(method3,wavelet3,level3,parameters,Ydata,ErrorType,zeta);
        parametershat_WW3e=fmincon(f,parameters0,[],[],[],[],LB,UB,[],options);
        
        % Store estimates in vectors
        d_vec_WW3(i)=parametershat_WW3e(1);
        omega_vec_WW3(i)=omg;
        a2_vec_WW3(i)=parametershat_WW3e(3);       
        b1_vec_WW3(i)=parametershat_WW3e(4);              
        theta_vec_WW3(i)=parametershat_WW3e(5);   
        delta_vec_WW3(i)=parametershat_WW3e(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;
   
        [~,eht_f,~]=FIEGARCH_core(Ydata_plus,d_vec_WW3(i),omega_vec_WW3(i),a2_vec_WW3(i),b1_vec_WW3(i),theta_vec_WW3(i),delta_vec_WW3(i), ErrorType,zeta);

        est_err_WW3=eht_f-ehtSim_plus;
        
        est_err_WW3_in_Q=est_err_WW3(end-insample-N_next+1:end-N_next);
        
        detectnansWW3=sum(isnan(est_err_WW3_in_Q))
        if detectnansWW3>0
         est_err_WW3_in_099=est_err_WW3_in_Q
        else
         est_err_WW3_in_099=est_err_WW3_in_Q(abs(est_err_WW3_in_Q)<quantile(abs(est_err_WW3_in_Q),0.99));
        end
        
        
        est_err_WW3_in(i)=mean(est_err_WW3_in_099);
        est_mad_WW3_in(i)=mean(abs(est_err_WW3_in_099));
        est_mse_WW3_in(i)=mean((est_err_WW3_in_099).^2);
        
        est_err_WW3_out(i)=mean(est_err_WW3(end-N_next+1:end));
        est_mad_WW3_out(i)=mean(abs(est_err_WW3(end-N_next+1:end)));
        est_mse_WW3_out(i)=mean((est_err_WW3(end-N_next+1:end)).^2);
        
        %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 % needs to be checked..before using
    
              JumpsSim='No Jumps';
              %JumpsSim='Jumps (0.028, N(0,0.2))';
              
              JumpsEst='Jumps Not Corrected';                            
              %JumpsEst='Jumps Corrected';
    
    disp('You are estimating zeta ')
    disp('Forecasting not yet available! ')
    %Extension not used in the thesis
 
    %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 
        % Detection in estimators! Turn on or off.. but
        % never forget to set what you really want
             
             [Ydata_plus,HtSim_plus,~]=FIEGARCH_simulate(N+N_next,d,omega,aL,bL,theta,delta,ErrorType,zeta); % Nx1
              Ydata=Ydata_plus(1:N);
              ehtSim_plus=exp(HtSim_plus);
              NextVals(i)=mean(ehtSim_plus(N+1:end))
              
          %{    
        % B) Jumps 
        % Detection in estimators! Turn on or off.. but
        % never forget to set what you really want
        
        % Adjust jump settings in the function
        
             [Ydata_plus,HtSim_plus,~,~]=FIEGARCH_Jump_simulate(N+N_next,d,omega,aL,bL,theta,delta,ErrorType,zeta); % Nx1
              Ydata=Ydata_plus(1:N);
              ehtSim_plus=exp(HtSim_plus);
              NextVals(i)=mean(ehtSim_plus(N+1:end));
        
              %}
              
        % omega estimate for FWE and WWE separately from other estimates 
         omg=mean(log(Ydata.^2))+1.2686;  % ln(yt^2)=ln(zt^2)+omega+inf mean zero MA
                                
                                 
        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_MLEe=fmincon(f_MLE,parameters0,[],[],[],[],LBMLE,UBMLE,[],options);
        
        % Store estimates in vectors
        d_vec_MLE(i)=parametershat_MLEe(1);
        omega_vec_MLE(i)=parametershat_MLEe(2);
        a2_vec_MLE(i)=parametershat_MLEe(3);       
        b1_vec_MLE(i)=parametershat_MLEe(4);              
        theta_vec_MLE(i)=parametershat_MLEe(5);   
        delta_vec_MLE(i)=parametershat_MLEe(6);  
            nu_vec_MLE(i)=NaN;
            skew_vec_MLE(i)=NaN;
        % Generate vectors for fixed parameters
        if strcmp(ErrorType,'stdged')
            nu_vec_MLE(i)=parametershat_MLEe(7);
            skew_vec_MLE(i)=NaN;
        elseif strcmp(ErrorType,'stdstudent')
            nu_vec_MLE(i)=parametershat_MLEe(7);
        skew_vec_MLE(i)=NaN;
        elseif strcmp(ErrorType,'skewstudent')
            nu_vec_MLE(i)=parametershat_MLEe(7);
            skew_vec_MLE(i)=parametershat_MLEe(8);
        end;
   
        
        % Estimate the eht series (conditional variance (not log-variance)
        %of the process) and calculate mean error for given simulation, mean absolute deviation and
        %mean squared error. This will be summed up later on and the overall RMSE will
        %be calculated (choose forecast interval carefully :) )
        
        [~,eht_f,~]=FIEGARCH_core(Ydata_plus,d_vec_MLE(i),omega_vec_MLE(i),a2_vec_MLE(i),b1_vec_MLE(i),theta_vec_MLE(i),delta_vec_MLE(i), ErrorType,[nu_vec_MLE(i),skew_vec_MLE(i)=parametershat_MLEe(8);]);
        
        est_err_MLE=eht_f-ehtSim_plus;
        
        est_err_MLE_in_Q=est_err_MLE(end-insample-N_next+1:end-N_next);
        
        detectnansMLE=sum(isnan(est_err_MLE_in_Q))
        if detectnansMLE>0
         est_err_MLE_in_099=est_err_MLE_in_Q
        else
         est_err_MLE_in_099=est_err_MLE_in_Q(abs(est_err_MLE_in_Q)<quantile(abs(est_err_MLE_in_Q),0.99));
        end
        
        est_err_MLE_in(i)=mean(est_err_MLE_in_099);
        est_mad_MLE_in(i)=mean(abs(est_err_MLE_in_099));
        est_mse_MLE_in(i)=mean((est_err_MLE_in_099).^2);
        
        est_err_MLE_out(i)=mean(est_err_MLE(end-N_next+1:end));
        est_mad_MLE_out(i)=mean(abs(est_err_MLE(end-N_next+1:end)));
        est_mse_MLE_out(i)=mean((est_err_MLE(end-N_next+1:end)).^2);
        
        %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_We=fmincon(f_W,parameters0,[],[],[],[],LB,UB,[],options);
        
        % Store estimates in vectors
        d_vec_W(i)=parametershat_We(1);
        omega_vec_W(i)=omg;
        a2_vec_W(i)=parametershat_We(3);       
        b1_vec_W(i)=parametershat_We(4);              
        theta_vec_W(i)=parametershat_We(5);   
        delta_vec_W(i)=parametershat_We(6);  
        nu_vec_W(i)=NaN;
        skew_vec_W(i)=NaN;
        
        % Generate vectors for fixed parameters
        if strcmp(ErrorType,'stdged')
            nu_vec_W(i)=parametershat_We(7); 
           skew_vec_W(i)=NaN;
        elseif strcmp(ErrorType,'stdstudent')
            nu_vec_W(i)=parametershat_We(7);
            skew_vec_W(i)=NaN;
        elseif strcmp(ErrorType,'skewstudent')
            nu_vec_W(i)=parametershat_We(7);
            skew_vec_W(i)=parametershat_We(8);
        end;
       
        [~,eht_f,~]=FIEGARCH_core(Ydata_plus,d_vec_W(i),omega_vec_W(i),a2_vec_W(i),b1_vec_W(i),theta_vec_W(i),delta_vec_W(i), ErrorType,[nu_vec_W(i),skew_vec_W(i)]);

        est_err_W=eht_f-ehtSim_plus;
        
        est_err_W_in_Q=est_err_W(end-insample-N_next+1:end-N_next);
        
        detectnansW=sum(isnan(est_err_W_in_Q))
        if detectnansW>0
         est_err_W_in_099=est_err_W_in_Q
        else
         est_err_W_in_099=est_err_W_in_Q(abs(est_err_W_in_Q)<quantile(abs(est_err_W_in_Q),0.99));
        end
        
   
        
        est_err_W_in(i)=mean(est_err_W_in_099);
        est_mad_W_in(i)=mean(abs(est_err_W_in_099));
        est_mse_W_in(i)=mean((est_err_W_in_099).^2);
        
        est_err_W_out(i)=mean(est_err_W(end-N_next+1:end));
        est_mad_W_out(i)=mean(abs(est_err_W(end-N_next+1:end)));
        est_mse_W_out(i)=mean((est_err_W(end-N_next+1:end)).^2);
        
        %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_WW1e=fmincon(f,parameters0,[],[],[],[],LB,UB,[],options);
        
        % Store estimates in vectors
        d_vec_WW1(i)=parametershat_WW1e(1);
        omega_vec_WW1(i)=omg;
        a2_vec_WW1(i)=parametershat_WW1e(3);       
        b1_vec_WW1(i)=parametershat_WW1e(4);              
        theta_vec_WW1(i)=parametershat_WW1e(5);   
        delta_vec_WW1(i)=parametershat_WW1e(6);  
        nu_vec_WW1(i)=NaN;
        skew_vec_WW1(i)=NaN;
        
        % Generate vectors for fixed parameters
        if strcmp(ErrorType,'stdged')
            nu_vec_WW1(i)=parametershat_WW1e(7);
            skew_vec_WW1(i)=NaN;
        elseif strcmp(ErrorType,'stdstudent')
            nu_vec_WW1(i)=parametershat_WW1e(7);
            skew_vec_WW1(i)=NaN;
        elseif strcmp(ErrorType,'skewstudent')
            nu_vec_WW1(i)=parametershat_WW1e(7);
            skew_vec_WW1(i)=parametershat_WW1e(8);
        end;
       
        [~,eht_f,~]=FIEGARCH_core(Ydata_plus,d_vec_WW1(i),omega_vec_WW1(i),a2_vec_WW1(i),b1_vec_WW1(i),theta_vec_WW1(i),delta_vec_WW1(i), ErrorType,[nu_vec_WW1(i),skew_vec_WW1(i)]);

       est_err_WW1=eht_f-ehtSim_plus;
        
        est_err_WW1_in_Q=est_err_WW1(end-insample-N_next+1:end-N_next);
        
        detectnansWW1=sum(isnan(est_err_WW1_in_Q))
        if detectnansWW1>0
         est_err_WW1_in_099=est_err_WW1_in_Q
        else
         est_err_WW1_in_099=est_err_WW1_in_Q(abs(est_err_WW1_in_Q)<quantile(abs(est_err_WW1_in_Q),0.99));
        end
        
        
        est_err_WW1_in(i)=mean(est_err_WW1_in_099);
        est_mad_WW1_in(i)=mean(abs(est_err_WW1_in_099));
        est_mse_WW1_in(i)=mean((est_err_WW1_in_099).^2);
        
        est_err_WW1_out(i)=mean(est_err_WW1(end-N_next+1:end));
        est_mad_WW1_out(i)=mean(abs(est_err_WW1(end-N_next+1:end)));
        est_mse_WW1_out(i)=mean((est_err_WW1(end-N_next+1:end)).^2);
        
        %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_WW2e=fmincon(f,parameters0,[],[],[],[],LB,UB,[],options);
        
        % Store estimates in vectors
        d_vec_WW2(i)=parametershat_WW2e(1);
        omega_vec_WW2(i)=omg;
        a2_vec_WW2(i)=parametershat_WW2e(3);       
        b1_vec_WW2(i)=parametershat_WW2e(4);              
        theta_vec_WW2(i)=parametershat_WW2e(5);   
        delta_vec_WW2(i)=parametershat_WW2e(6);  
        nu_vec_WW2(i)=NaN;
        skew_vec_WW2(i)=NaN;
        % Generate vectors for fixed parameters
        if strcmp(ErrorType,'stdged')
            nu_vec_WW2(i)=parametershat_WW2e(7);
            skew_vec_WW2(i)=NaN;
        elseif strcmp(ErrorType,'stdstudent')
            nu_vec_WW2(i)=parametershat_WW2e(7);
            skew_vec_WW2(i)=NaN;
        elseif strcmp(ErrorType,'skewstudent')
            nu_vec_WW2(i)=parametershat_WW2e(7);
            skew_vec_WW2(i)=parametershat_WW2e(8);
        end;

        [~,eht_f,~]=FIEGARCH_core(Ydata_plus,d_vec_WW2(i),omega_vec_WW2(i),a2_vec_WW2(i),b1_vec_WW2(i),theta_vec_WW2(i),delta_vec_WW2(i), ErrorType,[nu_vec_WW2(i),skew_vec_WW2(i)]);
         est_err_WW2=eht_f-ehtSim_plus;
        
        est_err_WW2_in_Q=est_err_WW2(end-insample-N_next+1:end-N_next);
        
         detectnansWW2=sum(isnan(est_err_WW2_in_Q))
        if detectnansWW2>0
         est_err_WW2_in_099=est_err_WW2_in_Q
        else
         est_err_WW2_in_099=est_err_WW2_in_Q(abs(est_err_WW2_in_Q)<quantile(abs(est_err_WW2_in_Q),0.99));
        end
        
        
        est_err_WW2_in(i)=mean(est_err_WW2_in_099);
        est_mad_WW2_in(i)=mean(abs(est_err_WW2_in_099));
        est_mse_WW2_in(i)=mean((est_err_WW2_in_099).^2);
        
        est_err_WW2_out(i)=mean(est_err_WW2(end-N_next+1:end));
        est_mad_WW2_out(i)=mean(abs(est_err_WW2(end-N_next+1:end)));
        est_mse_WW2_out(i)=mean((est_err_WW2(end-N_next+1:end)).^2);
  
        %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_WW3e=fmincon(f,parameters0,[],[],[],[],LB,UB,[],options);
        
        % Store estimates in vectors
        d_vec_WW3(i)=parametershat_WW3e(1);
        omega_vec_WW3(i)=omg;
        a2_vec_WW3(i)=parametershat_WW3e(3);       
        b1_vec_WW3(i)=parametershat_WW3e(4);              
        theta_vec_WW3(i)=parametershat_WW3e(5);   
        delta_vec_WW3(i)=parametershat_WW3e(6);
        nu_vec_WW3(i)=NaN;
        skew_vec_WW3(i)=NaN;
        
        % Generate vectors for fixed parameters
        if strcmp(ErrorType,'stdged')
            nu_vec_WW3(i)=parametershat_WW3e(7);
            skew_vec_WW3(i)=NaN;
        elseif strcmp(ErrorType,'stdstudent')
            nu_vec_WW3(i)=parametershat_WW3e(7);
            skew_vec_WW3(i)=NaN;
        elseif strcmp(ErrorType,'skewstudent')
            nu_vec_WW3(i)=parametershat_WW3e(7);
            skew_vec_WW3(i)=parametershat_WW3e(8);
        end;
    
        [~,eht_f,~]=FIEGARCH_core(Ydata_plus,d_vec_WW3(i),omega_vec_WW3(i),a2_vec_WW3(i),b1_vec_WW3(i),theta_vec_WW3(i),delta_vec_WW3(i), ErrorType,[nu_vec_WW3(i),skew_vec_WW3(i)]);
        est_err_WW3=eht_f-ehtSim_plus;
        
        est_err_WW3_in_Q=est_err_WW3(end-insample-N_next+1:end-N_next);
        
        detectnansWW3=sum(isnan(est_err_WW3_in_Q))
        if detectnansWW3>0
         est_err_WW3_in_099=est_err_WW3_in_Q
        else
         est_err_WW3_in_099=est_err_WW3_in_Q(abs(est_err_WW3_in_Q)<quantile(abs(est_err_WW3_in_Q),0.99));
        end
        
        
        est_err_WW3_in(i)=mean(est_err_WW3_in_099);
        est_mad_WW3_in(i)=mean(abs(est_err_WW3_in_099));
        est_mse_WW3_in(i)=mean((est_err_WW3_in_099).^2);
        
        est_err_WW3_out(i)=mean(est_err_WW3(end-N_next+1:end));
        est_mad_WW3_out(i)=mean(abs(est_err_WW3(end-N_next+1:end)));
        est_mse_WW3_out(i)=mean((est_err_WW3(end-N_next+1:end)).^2);
    
        
        %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

% Detect cases with poor convergence in fmincon and estim. at LB/UB, in
% general ommit estimates leading to huge errors or even NaNs in forecasted
% time series (+ count the number of cases when the forecasts are extremely wrong)

detectthres=(mean(NextVals)*10)*ones(M,1);% rule of thumb detection of extreme mean errors
% mle
detect_MLE1=abs(est_err_MLE_in)>detectthres;
detect_MLE2=isnan(est_err_MLE_in); % detect NaNs 
detect_MLE=(detect_MLE1+detect_MLE2)>0;% over threshold or NaNs?

valid_MLE=detect_MLE==0; % not over threshold, no NaNs?
fails_detected_MLE=sum(detect_MLE);% total 

N_NaN_MLE=sum(detect_MLE2);% num of NaNs

    % insample statistics for valid cases only
    est_err_MLE_in1=est_err_MLE_in(valid_MLE==1);
    est_mse_MLE_in1=est_mse_MLE_in(valid_MLE==1);
    est_mad_MLE_in1=est_mad_MLE_in(valid_MLE==1);
   
    est_err_MLE_out(valid_MLE==0)=inf;
    est_mse_MLE_out(valid_MLE==0)=inf;
    est_mad_MLE_out(valid_MLE==0)=inf; 
   
    est_err_MLE_out1= est_err_MLE_out(valid_MLE==1);
    est_mse_MLE_out1=est_mse_MLE_out(valid_MLE==1);
    est_mad_MLE_out1=est_mad_MLE_out(valid_MLE==1);
% w
detect_W1=abs(est_err_W_in)>detectthres;
detect_W2=isnan(est_err_W_in);
detect_W=(detect_W1+detect_W2)>0;

valid_W=detect_W==0;
fails_detected_W=sum(detect_W);

N_NaN_W=sum(detect_W2);

est_err_W_in1=est_err_W_in(valid_W==1);
est_mse_W_in1=est_mse_W_in(valid_W==1);
est_mad_W_in1=est_mad_W_in(valid_W==1);


    est_err_W_out(valid_W==0)=Inf;
    est_mse_W_out(valid_W==0)=Inf;
    est_mad_W_out(valid_W==0)=Inf ;

    est_err_W_out1= est_err_W_out(valid_W==1);
    est_mse_W_out1=est_mse_W_out(valid_W==1);
    est_mad_W_out1=est_mad_W_out(valid_W==1);


% ww1
detect_WW11=abs(est_err_WW1_in)>detectthres;
detect_WW12=isnan(est_err_WW1_in);
detect_WW1=(detect_WW11+detect_WW12)>0;

valid_WW1=detect_WW1==0;
fails_detected_WW1=sum(detect_WW1);

N_NaN_WW1=sum(detect_WW12);

est_err_WW1_in1=est_err_WW1_in(valid_WW1==1);
est_mse_WW1_in1=est_mse_WW1_in(valid_WW1==1);
est_mad_WW1_in1=est_mad_WW1_in(valid_WW1==1);


    est_err_WW1_out(valid_WW1==0)=inf;
    est_mse_WW1_out(valid_WW1==0)=inf;
    est_mad_WW1_out(valid_WW1==0)=inf; 

    est_err_WW1_out1= est_err_WW1_out(valid_WW1==1);
    est_mse_WW1_out1=est_mse_WW1_out(valid_WW1==1);
    est_mad_WW1_out1=est_mad_WW1_out(valid_WW1==1);

% ww2
detect_WW21=abs(est_err_WW2_in)>detectthres;
detect_WW22=isnan(est_err_WW2_in);
detect_WW2=(detect_WW21+detect_WW22)>0;

valid_WW2=detect_WW2==0;
fails_detected_WW2=sum(detect_WW2);

N_NaN_WW2=sum(detect_WW22);

est_err_WW2_in1=est_err_WW2_in(valid_WW2==1);
est_mse_WW2_in1=est_mse_WW2_in(valid_WW2==1);
est_mad_WW2_in1=est_mad_WW2_in(valid_WW2==1);


    est_err_WW2_out(valid_WW2==0)=inf;
    est_mse_WW2_out(valid_WW2==0)=inf;
    est_mad_WW2_out(valid_WW2==0)=inf; 

    est_err_WW2_out1= est_err_WW2_out(valid_WW2==1);
    est_mse_WW2_out1=est_mse_WW2_out(valid_WW2==1);
    est_mad_WW2_out1=est_mad_WW2_out(valid_WW2==1);

% ww3
detect_WW31=abs(est_err_WW3_in)>detectthres;
detect_WW32=isnan(est_err_WW3_in);
detect_WW3=(detect_WW31+detect_WW32)>0;

valid_WW3=detect_WW3==0;
fails_detected_WW3=sum(detect_WW3);

N_NaN_WW3=sum(detect_WW32);

est_err_WW3_in1=est_err_WW3_in(valid_WW3==1);
est_mse_WW3_in1=est_mse_WW3_in(valid_WW3==1);
est_mad_WW3_in1=est_mad_WW3_in(valid_WW3==1);


    est_err_WW3_out(valid_WW3==0)=inf;
    est_mse_WW3_out(valid_WW3==0)=inf;
    est_mad_WW3_out(valid_WW3==0)=inf; 
  
    est_err_WW3_out1= est_err_WW3_out(valid_WW3==1);
    est_mse_WW3_out1=est_mse_WW3_out(valid_WW3==1);
    est_mad_WW3_out1=est_mad_WW3_out(valid_WW3==1);
    
valid_all=((detect_MLE==1)+(detect_W==1)+(detect_WW1==1)+(detect_WW2==1)+(detect_WW3==1))==0;
N_valid_all=sum(valid_all);

est_err_MLE_in0=est_err_MLE_in(valid_all==1);
est_mse_MLE_in0=est_mse_MLE_in(valid_all==1);
est_mad_MLE_in0=est_mad_MLE_in(valid_all==1);
    est_err_MLE_out0= est_err_MLE_out(valid_all==1);
    est_mse_MLE_out0=est_mse_MLE_out(valid_all==1);
    est_mad_MLE_out0=est_mad_MLE_out(valid_all==1);
    
est_err_W_in0=est_err_W_in(valid_all==1);
est_mse_W_in0=est_mse_W_in(valid_all==1);
est_mad_W_in0=est_mad_W_in(valid_all==1);
    est_err_W_out0= est_err_W_out(valid_all==1);
    est_mse_W_out0=est_mse_W_out(valid_all==1);
    est_mad_W_out0=est_mad_W_out(valid_all==1);

est_err_WW1_in0=est_err_WW1_in(valid_all==1);
est_mse_WW1_in0=est_mse_WW1_in(valid_all==1);
est_mad_WW1_in0=est_mad_WW1_in(valid_all==1); 
    est_err_WW1_out0= est_err_WW1_out(valid_all==1);
    est_mse_WW1_out0=est_mse_WW1_out(valid_all==1);
    est_mad_WW1_out0=est_mad_WW1_out(valid_all==1);

est_err_WW2_in0=est_err_WW2_in(valid_all==1);
est_mse_WW2_in0=est_mse_WW2_in(valid_all==1);
est_mad_WW2_in0=est_mad_WW2_in(valid_all==1);
    est_err_WW2_out0= est_err_WW2_out(valid_all==1);
    est_mse_WW2_out0=est_mse_WW2_out(valid_all==1);
    est_mad_WW2_out0=est_mad_WW2_out(valid_all==1);

est_err_WW3_in0=est_err_WW3_in(valid_all==1);
est_mse_WW3_in0=est_mse_WW3_in(valid_all==1);
est_mad_WW3_in0=est_mad_WW3_in(valid_all==1);
    est_err_WW3_out0= est_err_WW3_out(valid_all==1);
    est_mse_WW3_out0=est_mse_WW3_out(valid_all==1);
    est_mad_WW3_out0=est_mad_WW3_out(valid_all==1);

%Main Statistics
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%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)); 

est_err_in_mean_MLE=mean(est_err_MLE_in1);
est_mad_in_mean_MLE=mean(est_mad_MLE_in1);
RMSE9_in_MLE=sqrt(mean(est_mse_MLE_in1)); 
est_err_out_mean_MLE=mean(est_err_MLE_out1);
est_mad_out_mean_MLE=mean(est_mad_MLE_out1);
RMSE9_out_MLE=sqrt(mean(est_mse_MLE_out1)); 

% individually corrected
Q_099_mad_out_MLE=quantile(est_mad_MLE_out1,0.99);
Q_095_mad_out_MLE=quantile(est_mad_MLE_out1,0.95);
Q_090_mad_out_MLE=quantile(est_mad_MLE_out1,0.90);
Q_050_mad_out_MLE=quantile(est_mad_MLE_out1,0.50);

Q_099_err_out_MLE=quantile(est_err_MLE_out1,0.99);
Q_095_err_out_MLE=quantile(est_err_MLE_out1,0.95);
Q_090_err_out_MLE=quantile(est_err_MLE_out1,0.90);
Q_010_err_out_MLE=quantile(est_err_MLE_out1,0.10);
Q_005_err_out_MLE=quantile(est_err_MLE_out1,0.05);
Q_001_err_out_MLE=quantile(est_err_MLE_out1,0.01);

%total
Q_099_mad_out_MLE_tot=quantile(est_mad_MLE_out,0.99);
Q_095_mad_out_MLE_tot=quantile(est_mad_MLE_out,0.95);
Q_090_mad_out_MLE_tot=quantile(est_mad_MLE_out,0.90);
Q_085_mad_out_MLE_tot=quantile(est_mad_MLE_out,0.85);
Q_080_mad_out_MLE_tot=quantile(est_mad_MLE_out,0.80);
Q_075_mad_out_MLE_tot=quantile(est_mad_MLE_out,0.75);
Q_070_mad_out_MLE_tot=quantile(est_mad_MLE_out,0.70);
Q_050_mad_out_MLE_tot=quantile(est_mad_MLE_out,0.50);

% equal simulation set
est_err_in_mean_MLE0=mean(est_err_MLE_in0);
est_mad_in_mean_MLE0=mean(est_mad_MLE_in0);
RMSE9_in_MLE0=sqrt(mean(est_mse_MLE_in0)); 
est_err_out_mean_MLE0=mean(est_err_MLE_out0);
est_mad_out_mean_MLE0=mean(est_mad_MLE_out0);
RMSE9_out_MLE0=sqrt(mean(est_mse_MLE_out0)); 

Q_099_mad_out_MLE0=quantile(est_mad_MLE_out0,0.99);
Q_095_mad_out_MLE0=quantile(est_mad_MLE_out0,0.95);
Q_090_mad_out_MLE0=quantile(est_mad_MLE_out0,0.90);
Q_050_mad_out_MLE0=quantile(est_mad_MLE_out0,0.50);

Q_099_err_out_MLE0=quantile(est_err_MLE_out0,0.99);
Q_095_err_out_MLE0=quantile(est_err_MLE_out0,0.95);
Q_090_err_out_MLE0=quantile(est_err_MLE_out0,0.90);
Q_010_err_out_MLE0=quantile(est_err_MLE_out0,0.10);
Q_005_err_out_MLE0=quantile(est_err_MLE_out0,0.05);
Q_001_err_out_MLE0=quantile(est_err_MLE_out0,0.01);


%corrected measures  for estimates providing "reasonable"forecasts( no NaNs, inf, huge mean absolute error)
Vdhat_MLE=mean(d_vec_MLE(valid_all==1));
Vbias1_MLE=Vdhat_MLE-d;                    
VRMSE1_MLE=sqrt(sum(((d_vec_MLE-d).*valid_all).^2)/N_valid_all);  

Vomegahat_MLE=mean(omega_vec_MLE(valid_all==1));
Vbias2_MLE=Vomegahat_MLE-omega;                    
VRMSE2_MLE=sqrt(sum(((omega_vec_MLE-omega).*valid_all).^2)/N_valid_all);  

Va2hat_MLE=mean(a2_vec_MLE(valid_all==1));                
Vbias3_MLE=Va2hat_MLE-a2;                    
VRMSE3_MLE=sqrt(sum(((a2_vec_MLE-a2).*valid_all).^2)/N_valid_all);  

Vb1hat_MLE=mean(b1_vec_MLE(valid_all==1));                
Vbias4_MLE=Vb1hat_MLE-b1;                   
VRMSE4_MLE=sqrt(sum(((b1_vec_MLE-b1).*valid_all).^2)/N_valid_all);  

Vthetahat_MLE=mean(theta_vec_MLE(valid_all==1));               
Vbias5_MLE=Vthetahat_MLE-theta;                  
VRMSE5_MLE=sqrt(sum(((theta_vec_MLE-theta).*valid_all).^2)/N_valid_all); 

Vdeltahat_MLE=mean(delta_vec_MLE(valid_all==1));              
Vbias6_MLE=Vdeltahat_MLE-delta;                  
VRMSE6_MLE=sqrt(sum(((delta_vec_MLE-delta).*valid_all).^2)/N_valid_all); 


% 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));

est_err_in_mean_W=mean(est_err_W_in1);
est_mad_in_mean_W=mean(est_mad_W_in1);
RMSE9_in_W=sqrt(mean(est_mse_W_in1)); 
est_err_out_mean_W=mean(est_err_W_out1);
est_mad_out_mean_W=mean(est_mad_W_out1);
RMSE9_out_W=sqrt(mean(est_mse_W_out1)); 

Q_099_mad_out_W=quantile(est_mad_W_out1,0.99);
Q_095_mad_out_W=quantile(est_mad_W_out1,0.95);
Q_090_mad_out_W=quantile(est_mad_W_out1,0.90);
Q_050_mad_out_W=quantile(est_mad_W_out1,0.50);

Q_099_err_out_W=quantile(est_err_W_out1,0.99);
Q_095_err_out_W=quantile(est_err_W_out1,0.95);
Q_090_err_out_W=quantile(est_err_W_out1,0.90);
Q_010_err_out_W=quantile(est_err_W_out1,0.10);
Q_005_err_out_W=quantile(est_err_W_out1,0.05);
Q_001_err_out_W=quantile(est_err_W_out1,0.01);

%total
Q_099_mad_out_W_tot=quantile(est_mad_W_out,0.99);
Q_095_mad_out_W_tot=quantile(est_mad_W_out,0.95);
Q_090_mad_out_W_tot=quantile(est_mad_W_out,0.90);
Q_085_mad_out_W_tot=quantile(est_mad_W_out,0.85);
Q_080_mad_out_W_tot=quantile(est_mad_W_out,0.80);
Q_075_mad_out_W_tot=quantile(est_mad_W_out,0.75);
Q_070_mad_out_W_tot=quantile(est_mad_W_out,0.70);
Q_050_mad_out_W_tot=quantile(est_mad_W_out,0.50);

% equal simulation set
est_err_in_mean_W0=mean(est_err_W_in0);
est_mad_in_mean_W0=mean(est_mad_W_in0);
RMSE9_in_W0=sqrt(mean(est_mse_W_in0)); 
est_err_out_mean_W0=mean(est_err_W_out0);
est_mad_out_mean_W0=mean(est_mad_W_out0);
RMSE9_out_W0=sqrt(mean(est_mse_W_out0)); 

Q_099_mad_out_W0=quantile(est_mad_W_out0,0.99);
Q_095_mad_out_W0=quantile(est_mad_W_out0,0.95);
Q_090_mad_out_W0=quantile(est_mad_W_out0,0.90);
Q_050_mad_out_W0=quantile(est_mad_W_out0,0.50);

Q_099_err_out_W0=quantile(est_err_W_out0,0.99);
Q_095_err_out_W0=quantile(est_err_W_out0,0.95);
Q_090_err_out_W0=quantile(est_err_W_out0,0.90);
Q_010_err_out_W0=quantile(est_err_W_out0,0.10);
Q_005_err_out_W0=quantile(est_err_W_out0,0.05);
Q_001_err_out_W0=quantile(est_err_W_out0,0.01);

%corrected measures for param
Vdhat_W=mean(d_vec_W(valid_all==1));
Vbias1_W=Vdhat_W-d;                    
VRMSE1_W=sqrt(sum(((d_vec_W-d).*valid_all).^2)/N_valid_all);  

Vomegahat_W=mean(omega_vec_W(valid_all==1));
Vbias2_W=Vomegahat_W-omega;                    
VRMSE2_W=sqrt(sum(((omega_vec_W-omega).*valid_all).^2)/N_valid_all);  

Va2hat_W=mean(a2_vec_W(valid_all==1));                
Vbias3_W=Va2hat_W-a2;                    
VRMSE3_W=sqrt(sum(((a2_vec_W-a2).*valid_all).^2)/N_valid_all);  

Vb1hat_W=mean(b1_vec_W(valid_all==1));                
Vbias4_W=Vb1hat_W-b1;                   
VRMSE4_W=sqrt(sum(((b1_vec_W-b1).*valid_all).^2)/N_valid_all);  

Vthetahat_W=mean(theta_vec_W(valid_all==1));               
Vbias5_W=Vthetahat_W-theta;                  
VRMSE5_W=sqrt(sum(((theta_vec_W-theta).*valid_all).^2)/N_valid_all); 

Vdeltahat_W=mean(delta_vec_W(valid_all==1));              
Vbias6_W=Vdeltahat_W-delta;                  
VRMSE6_W=sqrt(sum(((delta_vec_W-delta).*valid_all).^2)/N_valid_all); 

%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)); 

est_err_in_mean_WW1=mean(est_err_WW1_in1);
est_mad_in_mean_WW1=mean(est_mad_WW1_in1);
RMSE9_in_WW1=sqrt(mean(est_mse_WW1_in1)); 
est_err_out_mean_WW1=mean(est_err_WW1_out1);
est_mad_out_mean_WW1=mean(est_mad_WW1_out1);
RMSE9_out_WW1=sqrt(mean(est_mse_WW1_out1)); 

Q_099_mad_out_WW1=quantile(est_mad_WW1_out1,0.99);
Q_095_mad_out_WW1=quantile(est_mad_WW1_out1,0.95);
Q_090_mad_out_WW1=quantile(est_mad_WW1_out1,0.90);
Q_050_mad_out_WW1=quantile(est_mad_WW1_out1,0.50);

Q_099_err_out_WW1=quantile(est_err_WW1_out1,0.99);
Q_095_err_out_WW1=quantile(est_err_WW1_out1,0.95);
Q_090_err_out_WW1=quantile(est_err_WW1_out1,0.90);
Q_010_err_out_WW1=quantile(est_err_WW1_out1,0.10);
Q_005_err_out_WW1=quantile(est_err_WW1_out1,0.05);
Q_001_err_out_WW1=quantile(est_err_WW1_out1,0.01);
%total
Q_099_mad_out_WW1_tot=quantile(est_mad_WW1_out,0.99);
Q_095_mad_out_WW1_tot=quantile(est_mad_WW1_out,0.95);
Q_090_mad_out_WW1_tot=quantile(est_mad_WW1_out,0.90);
Q_085_mad_out_WW1_tot=quantile(est_mad_WW1_out,0.85);
Q_080_mad_out_WW1_tot=quantile(est_mad_WW1_out,0.80);
Q_075_mad_out_WW1_tot=quantile(est_mad_WW1_out,0.75);
Q_070_mad_out_WW1_tot=quantile(est_mad_WW1_out,0.70);
Q_050_mad_out_WW1_tot=quantile(est_mad_WW1_out,0.50);

% equal simulation set
est_err_in_mean_WW10=mean(est_err_WW1_in0);
est_mad_in_mean_WW10=mean(est_mad_WW1_in0);
RMSE9_in_WW10=sqrt(mean(est_mse_WW1_in0)); 
est_err_out_mean_WW10=mean(est_err_WW1_out0);
est_mad_out_mean_WW10=mean(est_mad_WW1_out0);
RMSE9_out_WW10=sqrt(mean(est_mse_WW1_out0)); 

Q_099_mad_out_WW10=quantile(est_mad_WW1_out0,0.99);
Q_095_mad_out_WW10=quantile(est_mad_WW1_out0,0.95);
Q_090_mad_out_WW10=quantile(est_mad_WW1_out0,0.90);
Q_050_mad_out_WW10=quantile(est_mad_WW1_out0,0.50);

Q_099_err_out_WW10=quantile(est_err_WW1_out0,0.99);
Q_095_err_out_WW10=quantile(est_err_WW1_out0,0.95);
Q_090_err_out_WW10=quantile(est_err_WW1_out0,0.90);
Q_010_err_out_WW10=quantile(est_err_WW1_out0,0.10);
Q_005_err_out_WW10=quantile(est_err_WW1_out0,0.05);
Q_001_err_out_WW10=quantile(est_err_WW1_out0,0.01);

%corrected measures for param
Vdhat_WW1=mean(d_vec_WW1(valid_all==1));
Vbias1_WW1=Vdhat_WW1-d;                    
VRMSE1_WW1=sqrt(sum(((d_vec_WW1-d).*valid_all).^2)/N_valid_all);  

Vomegahat_WW1=mean(omega_vec_WW1(valid_all==1));
Vbias2_WW1=Vomegahat_WW1-omega;                    
VRMSE2_WW1=sqrt(sum(((omega_vec_WW1-omega).*valid_all).^2)/N_valid_all);  

Va2hat_WW1=mean(a2_vec_WW1(valid_all==1));                
Vbias3_WW1=Va2hat_WW1-a2;                    
VRMSE3_WW1=sqrt(sum(((a2_vec_WW1-a2).*valid_all).^2)/N_valid_all);  

Vb1hat_WW1=mean(b1_vec_WW1(valid_all==1));                
Vbias4_WW1=Vb1hat_WW1-b1;                   
VRMSE4_WW1=sqrt(sum(((b1_vec_WW1-b1).*valid_all).^2)/N_valid_all);  

Vthetahat_WW1=mean(theta_vec_WW1(valid_all==1));               
Vbias5_WW1=Vthetahat_WW1-theta;                  
VRMSE5_WW1=sqrt(sum(((theta_vec_WW1-theta).*valid_all).^2)/N_valid_all); 

Vdeltahat_WW1=mean(delta_vec_WW1(valid_all==1));              
Vbias6_WW1=Vdeltahat_WW1-delta;                  
VRMSE6_WW1=sqrt(sum(((delta_vec_WW1-delta).*valid_all).^2)/N_valid_all); 

%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)); 

est_err_in_mean_WW2=mean(est_err_WW2_in1);
est_mad_in_mean_WW2=mean(est_mad_WW2_in1);
RMSE9_in_WW2=sqrt(mean(est_mse_WW2_in1)); 
est_err_out_mean_WW2=mean(est_err_WW2_out1);
est_mad_out_mean_WW2=mean(est_mad_WW2_out1);
RMSE9_out_WW2=sqrt(mean(est_mse_WW2_out1)); 

Q_099_mad_out_WW2=quantile(est_mad_WW2_out1,0.99);
Q_095_mad_out_WW2=quantile(est_mad_WW2_out1,0.95);
Q_090_mad_out_WW2=quantile(est_mad_WW2_out1,0.90);
Q_050_mad_out_WW2=quantile(est_mad_WW2_out1,0.50);

Q_099_err_out_WW2=quantile(est_err_WW2_out1,0.99);
Q_095_err_out_WW2=quantile(est_err_WW2_out1,0.95);
Q_090_err_out_WW2=quantile(est_err_WW2_out1,0.90);
Q_010_err_out_WW2=quantile(est_err_WW2_out1,0.10);
Q_005_err_out_WW2=quantile(est_err_WW2_out1,0.05);
Q_001_err_out_WW2=quantile(est_err_WW2_out1,0.01);

%total
Q_099_mad_out_WW2_tot=quantile(est_mad_WW2_out,0.99);
Q_095_mad_out_WW2_tot=quantile(est_mad_WW2_out,0.95);
Q_090_mad_out_WW2_tot=quantile(est_mad_WW2_out,0.90);
Q_085_mad_out_WW2_tot=quantile(est_mad_WW2_out,0.85);
Q_080_mad_out_WW2_tot=quantile(est_mad_WW2_out,0.80);
Q_075_mad_out_WW2_tot=quantile(est_mad_WW2_out,0.75);
Q_070_mad_out_WW2_tot=quantile(est_mad_WW2_out,0.70);
Q_050_mad_out_WW2_tot=quantile(est_mad_WW2_out,0.50);

% equal simulation set
est_err_in_mean_WW20=mean(est_err_WW2_in0);
est_mad_in_mean_WW20=mean(est_mad_WW2_in0);
RMSE9_in_WW20=sqrt(mean(est_mse_WW2_in0)); 
est_err_out_mean_WW20=mean(est_err_WW2_out0);
est_mad_out_mean_WW20=mean(est_mad_WW2_out0);
RMSE9_out_WW20=sqrt(mean(est_mse_WW2_out0)); 

Q_099_mad_out_WW20=quantile(est_mad_WW2_out0,0.99);
Q_095_mad_out_WW20=quantile(est_mad_WW2_out0,0.95);
Q_090_mad_out_WW20=quantile(est_mad_WW2_out0,0.90);
Q_050_mad_out_WW20=quantile(est_mad_WW2_out0,0.50);

Q_099_err_out_WW20=quantile(est_err_WW2_out0,0.99);
Q_095_err_out_WW20=quantile(est_err_WW2_out0,0.95);
Q_090_err_out_WW20=quantile(est_err_WW2_out0,0.90);
Q_010_err_out_WW20=quantile(est_err_WW2_out0,0.10);
Q_005_err_out_WW20=quantile(est_err_WW2_out0,0.05);
Q_001_err_out_WW20=quantile(est_err_WW2_out0,0.01);

%corrected measures for param
Vdhat_WW2=mean(d_vec_WW2(valid_all==1));
Vbias1_WW2=Vdhat_WW2-d;                    
VRMSE1_WW2=sqrt(sum(((d_vec_WW2-d).*valid_all).^2)/N_valid_all);  

Vomegahat_WW2=mean(omega_vec_WW2(valid_all==1));
Vbias2_WW2=Vomegahat_WW2-omega;                    
VRMSE2_WW2=sqrt(sum(((omega_vec_WW2-omega).*valid_all).^2)/N_valid_all);  

Va2hat_WW2=mean(a2_vec_WW2(valid_all==1));                
Vbias3_WW2=Va2hat_WW2-a2;                    
VRMSE3_WW2=sqrt(sum(((a2_vec_WW2-a2).*valid_all).^2)/N_valid_all);  

Vb1hat_WW2=mean(b1_vec_WW2(valid_all==1));                
Vbias4_WW2=Vb1hat_WW2-b1;                   
VRMSE4_WW2=sqrt(sum(((b1_vec_WW2-b1).*valid_all).^2)/N_valid_all);  

Vthetahat_WW2=mean(theta_vec_WW2(valid_all==1));               
Vbias5_WW2=Vthetahat_WW2-theta;                  
VRMSE5_WW2=sqrt(sum(((theta_vec_WW2-theta).*valid_all).^2)/N_valid_all); 

Vdeltahat_WW2=mean(delta_vec_WW2(valid_all==1));              
Vbias6_WW2=Vdeltahat_WW2-delta;                  
VRMSE6_WW2=sqrt(sum(((delta_vec_WW2-delta).*valid_all).^2)/N_valid_all); 


%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)); 

est_err_in_mean_WW3=mean(est_err_WW3_in1);
est_mad_in_mean_WW3=mean(est_mad_WW3_in1);
RMSE9_in_WW3=sqrt(mean(est_mse_WW3_in1)); 
est_err_out_mean_WW3=mean(est_err_WW3_out1);
est_mad_out_mean_WW3=mean(est_mad_WW3_out1);
RMSE9_out_WW3=sqrt(mean(est_mse_WW3_out1));

Q_099_mad_out_WW3=quantile(est_mad_WW3_out1,0.99);
Q_095_mad_out_WW3=quantile(est_mad_WW3_out1,0.95);
Q_090_mad_out_WW3=quantile(est_mad_WW3_out1,0.90);
Q_050_mad_out_WW3=quantile(est_mad_WW3_out1,0.50);

Q_099_err_out_WW3=quantile(est_err_WW3_out1,0.99);
Q_095_err_out_WW3=quantile(est_err_WW3_out1,0.95);
Q_090_err_out_WW3=quantile(est_err_WW3_out1,0.90);
Q_010_err_out_WW3=quantile(est_err_WW3_out1,0.10);
Q_005_err_out_WW3=quantile(est_err_WW3_out1,0.05);
Q_001_err_out_WW3=quantile(est_err_WW3_out1,0.01);

%total
Q_099_mad_out_WW3_tot=quantile(est_mad_WW3_out,0.99);
Q_095_mad_out_WW3_tot=quantile(est_mad_WW3_out,0.95);
Q_090_mad_out_WW3_tot=quantile(est_mad_WW3_out,0.90);
Q_085_mad_out_WW3_tot=quantile(est_mad_WW3_out,0.85);
Q_080_mad_out_WW3_tot=quantile(est_mad_WW3_out,0.80);
Q_075_mad_out_WW3_tot=quantile(est_mad_WW3_out,0.75);
Q_070_mad_out_WW3_tot=quantile(est_mad_WW3_out,0.70);
Q_050_mad_out_WW3_tot=quantile(est_mad_WW3_out,0.50);

% equal simulation set
est_err_in_mean_WW30=mean(est_err_WW3_in0);
est_mad_in_mean_WW30=mean(est_mad_WW3_in0);
RMSE9_in_WW30=sqrt(mean(est_mse_WW3_in0)); 
est_err_out_mean_WW30=mean(est_err_WW3_out0);
est_mad_out_mean_WW30=mean(est_mad_WW3_out0);
RMSE9_out_WW30=sqrt(mean(est_mse_WW3_out0)); 

Q_099_mad_out_WW30=quantile(est_mad_WW3_out0,0.99);
Q_095_mad_out_WW30=quantile(est_mad_WW3_out0,0.95);
Q_090_mad_out_WW30=quantile(est_mad_WW3_out0,0.90);
Q_050_mad_out_WW30=quantile(est_mad_WW3_out0,0.50);

Q_099_err_out_WW30=quantile(est_err_WW3_out0,0.99);
Q_095_err_out_WW30=quantile(est_err_WW3_out0,0.95);
Q_090_err_out_WW30=quantile(est_err_WW3_out0,0.90);
Q_010_err_out_WW30=quantile(est_err_WW3_out0,0.10);
Q_005_err_out_WW30=quantile(est_err_WW3_out0,0.05);
Q_001_err_out_WW30=quantile(est_err_WW3_out0,0.01);

%corrected measures for param
Vdhat_WW3=mean(d_vec_WW3(valid_all==1));
Vbias1_WW3=Vdhat_WW3-d;                    
VRMSE1_WW3=sqrt(sum(((d_vec_WW3-d).*valid_all).^2)/N_valid_all);  

Vomegahat_WW3=mean(omega_vec_WW3(valid_all==1));
Vbias2_WW3=Vomegahat_WW3-omega;                    
VRMSE2_WW3=sqrt(sum(((omega_vec_WW3-omega).*valid_all).^2)/N_valid_all);  

Va2hat_WW3=mean(a2_vec_WW3(valid_all==1));                
Vbias3_WW3=Va2hat_WW3-a2;                    
VRMSE3_WW3=sqrt(sum(((a2_vec_WW3-a2).*valid_all).^2)/N_valid_all);  

Vb1hat_WW3=mean(b1_vec_WW3(valid_all==1));                
Vbias4_WW3=Vb1hat_WW3-b1;                   
VRMSE4_WW3=sqrt(sum(((b1_vec_WW3-b1).*valid_all).^2)/N_valid_all);  

Vthetahat_WW3=mean(theta_vec_WW3(valid_all==1));               
Vbias5_WW3=Vthetahat_WW3-theta;                  
VRMSE5_WW3=sqrt(sum(((theta_vec_WW3-theta).*valid_all).^2)/N_valid_all); 

Vdeltahat_WW3=mean(delta_vec_WW3(valid_all==1));              
Vbias6_WW3=Vdeltahat_WW3-delta;                  
VRMSE6_WW3=sqrt(sum(((delta_vec_WW3-delta).*valid_all).^2)/N_valid_all); 


 %Summary
if strcmp(ErrorType,'stdnormal') % corrected values nonaffected by diverding cases (estimates at boundaries for almost all coefficients... )
    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=[];
    
    %corrected
    Vparametershat_MLE=[Vdhat_MLE Vomegahat_MLE Va2hat_MLE Vb1hat_MLE Vthetahat_MLE Vdeltahat_MLE  ];
    Vbias_MLE=[Vbias1_MLE Vbias2_MLE Vbias3_MLE Vbias4_MLE Vbias5_MLE Vbias6_MLE ];
    VRMSE_MLE=[VRMSE1_MLE VRMSE2_MLE VRMSE3_MLE VRMSE4_MLE VRMSE5_MLE VRMSE6_MLE ];
    
    
    %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=[];
    
    %corrected
    Vparametershat_W=[Vdhat_W Vomegahat_W Va2hat_W Vb1hat_W Vthetahat_W Vdeltahat_W ];
    Vbias_W=[Vbias1_W Vbias2_W Vbias3_W Vbias4_W Vbias5_W Vbias6_W ];
    VRMSE_W=[VRMSE1_W VRMSE2_W VRMSE3_W VRMSE4_W VRMSE5_W VRMSE6_W ];


    %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=[];
    
    %corrected
    Vparametershat_WW1=[Vdhat_WW1 Vomegahat_WW1 Va2hat_WW1 Vb1hat_WW1 Vthetahat_WW1 Vdeltahat_WW1 ];
    Vbias_WW1=[Vbias1_WW1 Vbias2_WW1 Vbias3_WW1 Vbias4_WW1 Vbias5_WW1 Vbias6_WW1 ];
    VRMSE_WW1=[VRMSE1_WW1 VRMSE2_WW1 VRMSE3_WW1 VRMSE4_WW1 VRMSE5_WW1 VRMSE6_WW1 ];

    
    %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];
    nu=[];
    
    %corrected
    Vparametershat_WW2=[Vdhat_WW2 Vomegahat_WW2 Va2hat_WW2 Vb1hat_WW2 Vthetahat_WW2 Vdeltahat_WW2];
    Vbias_WW2=[Vbias1_WW2 Vbias2_WW2 Vbias3_WW2 Vbias4_WW2 Vbias5_WW2 Vbias6_WW2 ];
    VRMSE_WW2=[VRMSE1_WW2 VRMSE2_WW2 VRMSE3_WW2 VRMSE4_WW2 VRMSE5_WW2 VRMSE6_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=[];
    
    %corrected
    Vparametershat_WW3=[Vdhat_WW3 Vomegahat_WW3 Va2hat_WW3 Vb1hat_WW3 Vthetahat_WW3 Vdeltahat_WW3 ];
    Vbias_WW3=[Vbias1_WW3 Vbias2_WW3 Vbias3_WW3 Vbias4_WW3 Vbias5_WW3 Vbias6_WW3 ];
    VRMSE_WW3=[VRMSE1_WW3 VRMSE2_WW3 VRMSE3_WW3 VRMSE4_WW3 VRMSE5_WW3 VRMSE6_WW3  ];
    
    
elseif strcmp(ErrorType,'stdged')
    Row_names={'\hat{d}}' '\hat{omega}' '\hat{a2}' '\hat{b1}' '\hat{theta}' '\hat{delta}' '\hat{nu}' 'forecast'};
    %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}' 'forecast'};
    %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 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}' 'forecast'};
     %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));
 disp(level1)
 disp(level2)
 disp(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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Corrected table:
diary off;

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)=Vparametershat_WW1(t);
    rmse(i)=VRMSE_WW1(t);
    bias(i)=Vbias_WW1(t);
    
    
    i=i+1;
   % param(i)='';
   
     hat(i)=Vparametershat_WW2(t);
    rmse(i)=VRMSE_WW2(t);
    bias(i)=Vbias_WW2(t);
    i=i+1;
    %param(i)='';
    
     hat(i)=Vparametershat_WW3(t);
    rmse(i)=VRMSE_WW3(t);
    bias(i)=Vbias_WW3(t);
    i=i+1;
    %param(i)='';
    
     hat(i)=Vparametershat_W(t);
    rmse(i)=VRMSE_W(t);
    bias(i)=Vbias_W(t);
    i=i+1;
   % param(i)='';
  
     hat(i)=Vparametershat_MLE(t);
    rmse(i)=VRMSE_MLE(t);
    bias(i)=Vbias_MLE(t);
end

diary Diary_file %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
diary on;
disp('Corrected table: Not influenced by not converging fmincon')
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('Corrected simulation set summary')
disp(['MLE insample:&mean error & mean abs.deviation & rmse &  =  &  ',num2str(est_err_in_mean_MLE),'  &  ',num2str(est_mad_in_mean_MLE),'  &  ',num2str(RMSE9_in_MLE),'  &']);
disp(['MLE outsample:&mean error & mean abs.deviation & rmse &  =  &  ',num2str(est_err_out_mean_MLE),'  &  ',num2str(est_mad_out_mean_MLE),'  &  ',num2str(RMSE9_out_MLE),'  &']);

disp(['W insample:&mean error & mean abs.deviation & rmse &  =  &  ',num2str(est_err_in_mean_W),'  &  ',num2str(est_mad_in_mean_W),'  &  ',num2str(RMSE9_in_W),'  &']);
disp(['W outsample:&mean error & mean abs.deviation & rmse &  =  &  ',num2str(est_err_out_mean_W),'  &  ',num2str(est_mad_out_mean_W),'  &  ',num2str(RMSE9_out_W),'  &']);

disp(['WW1 insample:&mean error & mean abs.deviation & rmse &  =  &  ',num2str(est_err_in_mean_WW1),'  &  ',num2str(est_mad_in_mean_WW1),'  &  ',num2str(RMSE9_in_WW1),'  &']);
disp(['WW1 outsample:&mean error & mean abs.deviation & rmse &  =  &  ',num2str(est_err_out_mean_WW1),'  &  ',num2str(est_mad_out_mean_WW1),'  &  ',num2str(RMSE9_out_WW1),'  &']);

disp(['WW2 insample:&mean error & mean abs.deviation & rmse &  =  &  ',num2str(est_err_in_mean_WW2),'  &  ',num2str(est_mad_in_mean_WW2),'  &  ',num2str(RMSE9_in_WW2),'  &']);
disp(['WW2 outsample:&mean error & mean abs.deviation & rmse &  =  &  ',num2str(est_err_out_mean_WW2),'  &  ',num2str(est_mad_out_mean_WW2),'  &  ',num2str(RMSE9_out_WW2),'  &']);

disp(['WW3 insample:&mean error & mean abs.deviation & rmse &  =  &  ',num2str(est_err_in_mean_WW3),'  &  ',num2str(est_mad_in_mean_WW3),'  &  ',num2str(RMSE9_in_WW3),'  &']);
disp(['WW3 outsample:&mean error & mean abs.deviation & rmse &  =  &  ',num2str(est_err_out_mean_WW3),'  &  ',num2str(est_mad_out_mean_WW3),'  &  ',num2str(RMSE9_out_WW3),'  &']);

% Outsample MAD Quantiles %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
disp('Corrected simulation set out Mad')
disp(['MLE mad outsample quantiles:& 0.50 & 0.90 & 0.95 & 0.99 &  =  &  ',num2str(Q_050_mad_out_MLE),'  &',num2str(Q_090_mad_out_MLE),'  &  ',num2str(Q_095_mad_out_MLE),'  &  ',num2str(Q_099_mad_out_MLE),'  &']);

disp(['W  mad outsample quantiles:& 0.50 & 0.90 &0.95 & 0.99 &  =  & ',num2str(Q_050_mad_out_W),'  & ',num2str(Q_090_mad_out_W),'  &  ',num2str(Q_095_mad_out_W),'  &  ',num2str(Q_099_mad_out_W),'  &']);

disp(['WW1  mad outsample quantiles:& 0.50 & 0.90 &0.95 & 0.99 &  = & ',num2str(Q_050_mad_out_WW1),'  &  ',num2str(Q_090_mad_out_WW1),'  &  ',num2str(Q_095_mad_out_WW1),'  &  ',num2str(Q_099_mad_out_WW1),'  &']);

disp(['WW2  mad outsample quantiles:& 0.50 &0.90 &0.95 & 0.99 &  =  & ',num2str(Q_050_mad_out_WW2),'  &  ',num2str(Q_090_mad_out_WW2),'  &  ',num2str(Q_095_mad_out_WW2),'  &  ',num2str(Q_099_mad_out_WW2),'  &']);

disp(['WW3  mad outsample quantiles:& 0.50 &0.90 &0.95 & 0.99 &  =  & ',num2str(Q_050_mad_out_WW3),'  & ',num2str(Q_090_mad_out_WW3),'  &  ',num2str(Q_095_mad_out_WW3),'  &  ',num2str(Q_099_mad_out_WW3),'  &']);

% Outsample RMSE Quantiles: currently N_next=1, thus
% for ith observation, i=1,...M:
% sqrt(sum(est_-true).^2/N_next)=sqrt((est-true).^2/1)=abs(est-true)=mad=|est-true|/1

% Outsample Error Quantiles Corrected %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp('Corrected simulation set out Err')
disp(['MLE err outsample quantiles:&0.90 &0.95 & 0.99 & =  &  ',num2str(Q_090_err_out_MLE),'  &  ',num2str(Q_095_err_out_MLE),'  &  ',num2str(Q_099_err_out_MLE),'  &']);
disp(['MLE err outsample quantiles:&0.01 &0.05 & 0.10 & =  &  ',num2str(Q_001_err_out_MLE),'  &  ',num2str(Q_005_err_out_MLE),'  &  ',num2str(Q_010_err_out_MLE),'  &']);

disp(['W err outsample quantiles:&0.90 &0.95 & 0.99 &  =  &  ',num2str(Q_090_err_out_W),'  &  ',num2str(Q_095_err_out_W),'  &  ',num2str(Q_099_err_out_W),'  &']);
disp(['W err outsample quantiles:&0.01 &0.05 & 0.10 &  =  &  ',num2str(Q_001_err_out_W),'  &  ',num2str(Q_005_err_out_W),'  &  ',num2str(Q_010_err_out_W),'  &']);

disp(['WW1 err outsample quantiles:&0.90 &0.95 & 0.99 &  =  &  ',num2str(Q_090_err_out_WW1),'  &  ',num2str(Q_095_err_out_WW1),'  &  ',num2str(Q_099_err_out_WW1),'  &']);
disp(['WW1 err outsample quantiles:&0.01 &0.05 & 0.10 &  =  &  ',num2str(Q_001_err_out_WW1),'  &  ',num2str(Q_005_err_out_WW1),'  &  ',num2str(Q_010_err_out_WW1),'  &']);

disp(['WW2 err outsample quantiles:&0.90 &0.95 & 0.99 &  =  &  ',num2str(Q_090_err_out_WW2),'  &  ',num2str(Q_095_err_out_WW2),'  &  ',num2str(Q_099_err_out_WW2),'  &']);
disp(['WW2 err outsample quantiles:&0.01 &0.05 & 0.10 &  =  &  ',num2str(Q_001_err_out_WW2),'  &  ',num2str(Q_005_err_out_WW2),'  &  ',num2str(Q_010_err_out_WW2),'  &']);

disp(['WW3 err outsample quantiles:&0.90 &0.95 & 0.99 &  =  &  ',num2str(Q_090_err_out_WW3),'  &  ',num2str(Q_095_err_out_WW3),'  &  ',num2str(Q_099_err_out_WW3),'  &']);
disp(['WW3 err outsample quantiles:&0.01 &0.05 & 0.10 &  =  &  ',num2str(Q_001_err_out_WW3),'  &  ',num2str(Q_005_err_out_WW3),'  &  ',num2str(Q_010_err_out_WW3),'  &']);

disp('End Corrected simulation set')
% Outsample MAD Quantiles Total %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
disp('Total simulation set out Mad')
disp(['MLE mad outsample quantiles total:& 0.50 & 0.70 & 0.75 & 0.80 & 0.85 & 0.90 & 0.95 & 0.99 &  =  &  ',num2str(Q_050_mad_out_MLE_tot),'  &',num2str(Q_070_mad_out_MLE_tot),'  &  ',num2str(Q_075_mad_out_MLE_tot),'  &',num2str(Q_080_mad_out_MLE_tot),'  &  ',num2str(Q_085_mad_out_MLE_tot),'  &',num2str(Q_090_mad_out_MLE_tot),'  &  ',num2str(Q_095_mad_out_MLE_tot),'  &  ',num2str(Q_099_mad_out_MLE_tot),'  &']);

disp(['W  mad outsample quantiles total:& 0.50 & 0.70 & 0.75 & 0.80 & 0.85 & 0.90 & 0.95 & 0.99 & =  & ',num2str(Q_050_mad_out_W_tot),'  & ',num2str(Q_070_mad_out_W_tot),'  &  ',num2str(Q_075_mad_out_W_tot),'  & ',num2str(Q_080_mad_out_W_tot),'  &  ',num2str(Q_085_mad_out_W_tot),'  & ',num2str(Q_090_mad_out_W_tot),'  &  ',num2str(Q_095_mad_out_W_tot),'  &  ',num2str(Q_099_mad_out_W_tot),'  &']);

disp(['WW1  mad outsample quantiles total:& 0.50 & 0.70 & 0.75 & 0.80 & 0.85 & 0.90 & 0.95 & 0.99 & = & ',num2str(Q_050_mad_out_WW1_tot),'  &  ',num2str(Q_070_mad_out_WW1_tot),'  &  ',num2str(Q_075_mad_out_WW1_tot),'  &  ',num2str(Q_080_mad_out_WW1_tot),'  &  ',num2str(Q_085_mad_out_WW1_tot),'   & ',num2str(Q_090_mad_out_WW1_tot),'  &  ',num2str(Q_095_mad_out_WW1_tot),'  &  ',num2str(Q_099_mad_out_WW1_tot),'  &']);

disp(['WW2  mad outsample quantiles total:& 0.50 & 0.70 & 0.75 & 0.80 & 0.85 & 0.90 & 0.95 & 0.99 &  =  & ',num2str(Q_050_mad_out_WW2_tot),'  &  ',num2str(Q_070_mad_out_WW2_tot),'  &  ',num2str(Q_075_mad_out_WW2_tot),'  &  ',num2str(Q_080_mad_out_WW2_tot),'  &  ',num2str(Q_085_mad_out_WW2_tot),'  &  ',num2str(Q_090_mad_out_WW2_tot),'  &  ',num2str(Q_095_mad_out_WW2_tot),'  &  ',num2str(Q_099_mad_out_WW2_tot),'  &']);

disp(['WW3  mad outsample quantiles total:& 0.50 & 0.70 & 0.75 & 0.80 & 0.85 & 0.90 & 0.95 & 0.99 &  =  & ',num2str(Q_050_mad_out_WW3_tot),'  & ',num2str(Q_070_mad_out_WW3_tot),'  &  ',num2str(Q_075_mad_out_WW3_tot),'  & ',num2str(Q_080_mad_out_WW3_tot),'  &  ',num2str(Q_085_mad_out_WW3_tot),'  & ',num2str(Q_090_mad_out_WW3_tot),'  &  ',num2str(Q_095_mad_out_WW3_tot),'  &  ',num2str(Q_099_mad_out_WW3_tot),'  &']);

disp('End Total simulation set ')

% Equal simulation set
disp('Equal simulation set summary')
disp(['MLE insample:&mean error & mean abs.deviation & rmse &  =  &  ',num2str(est_err_in_mean_MLE0),'  &  ',num2str(est_mad_in_mean_MLE0),'  &  ',num2str(RMSE9_in_MLE0),'  &']);
disp(['MLE outsample:&mean error & mean abs.deviation & rmse &  =  &  ',num2str(est_err_out_mean_MLE0),'  &  ',num2str(est_mad_out_mean_MLE0),'  &  ',num2str(RMSE9_out_MLE0),'  &']);

disp(['W insample:&mean error & mean abs.deviation & rmse &  =  &  ',num2str(est_err_in_mean_W0),'  &  ',num2str(est_mad_in_mean_W0),'  &  ',num2str(RMSE9_in_W0),'  &']);
disp(['W outsample:&mean error & mean abs.deviation & rmse &  =  &  ',num2str(est_err_out_mean_W0),'  &  ',num2str(est_mad_out_mean_W0),'  &  ',num2str(RMSE9_out_W0),'  &']);

disp(['WW1 insample:&mean error & mean abs.deviation & rmse &  =  &  ',num2str(est_err_in_mean_WW10),'  &  ',num2str(est_mad_in_mean_WW10),'  &  ',num2str(RMSE9_in_WW10),'  &']);
disp(['WW1 outsample:&mean error & mean abs.deviation & rmse &  =  &  ',num2str(est_err_out_mean_WW10),'  &  ',num2str(est_mad_out_mean_WW10),'  &  ',num2str(RMSE9_out_WW10),'  &']);

disp(['WW2 insample:&mean error & mean abs.deviation & rmse &  =  &  ',num2str(est_err_in_mean_WW20),'  &  ',num2str(est_mad_in_mean_WW20),'  &  ',num2str(RMSE9_in_WW20),'  &']);
disp(['WW2 outsample:&mean error & mean abs.deviation & rmse &  =  &  ',num2str(est_err_out_mean_WW20),'  &  ',num2str(est_mad_out_mean_WW20),'  &  ',num2str(RMSE9_out_WW20),'  &']);

disp(['WW3 insample:&mean error & mean abs.deviation & rmse &  =  &  ',num2str(est_err_in_mean_WW30),'  &  ',num2str(est_mad_in_mean_WW30),'  &  ',num2str(RMSE9_in_WW30),'  &']);
disp(['WW3 outsample:&mean error & mean abs.deviation & rmse &  =  &  ',num2str(est_err_out_mean_WW30),'  &  ',num2str(est_mad_out_mean_WW30),'  &  ',num2str(RMSE9_out_WW30),'  &']);

% Outsample MAD Quantiles %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
disp('Equal simulation set out Mad')
disp(['MLE mad outsample quantiles:& 0.50 & 0.90 & 0.95 & 0.99 &  =  &  ',num2str(Q_050_mad_out_MLE0),'  &',num2str(Q_090_mad_out_MLE0),'  &  ',num2str(Q_095_mad_out_MLE0),'  &  ',num2str(Q_099_mad_out_MLE0),'  &']);

disp(['W  mad outsample quantiles:& 0.50 & 0.90 &0.95 & 0.99 &  =  & ',num2str(Q_050_mad_out_W0),'  & ',num2str(Q_090_mad_out_W0),'  &  ',num2str(Q_095_mad_out_W0),'  &  ',num2str(Q_099_mad_out_W0),'  &']);

disp(['WW1  mad outsample quantiles:& 0.50 & 0.90 &0.95 & 0.99 &  = & ',num2str(Q_050_mad_out_WW10),'  &  ',num2str(Q_090_mad_out_WW10),'  &  ',num2str(Q_095_mad_out_WW10),'  &  ',num2str(Q_099_mad_out_WW10),'  &']);

disp(['WW2  mad outsample quantiles:& 0.50 &0.90 &0.95 & 0.99 &  =  & ',num2str(Q_050_mad_out_WW20),'  &  ',num2str(Q_090_mad_out_WW20),'  &  ',num2str(Q_095_mad_out_WW20),'  &  ',num2str(Q_099_mad_out_WW20),'  &']);

disp(['WW3  mad outsample quantiles:& 0.50 &0.90 &0.95 & 0.99 &  =  & ',num2str(Q_050_mad_out_WW30),'  & ',num2str(Q_090_mad_out_WW30),'  &  ',num2str(Q_095_mad_out_WW30),'  &  ',num2str(Q_099_mad_out_WW30),'  &']);

% Outsample RMSE Quantiles: currently N_next=1, thus
% for ith observation, i=1,...M:
% sqrt(sum(est_-true).^2/N_next)=sqrt((est-true).^2/1)=abs(est-true)=mad=|est-true|/1

% Outsample Error Quantiles Corrected %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp('Equal simulation set out Err')
disp(['MLE err outsample quantiles:&0.90 &0.95 & 0.99 & =  &  ',num2str(Q_090_err_out_MLE0),'  &  ',num2str(Q_095_err_out_MLE0),'  &  ',num2str(Q_099_err_out_MLE0),'  &']);
disp(['MLE err outsample quantiles:&0.01 &0.05 & 0.10 & =  &  ',num2str(Q_001_err_out_MLE0),'  &  ',num2str(Q_005_err_out_MLE0),'  &  ',num2str(Q_010_err_out_MLE0),'  &']);

disp(['W err outsample quantiles:&0.90 &0.95 & 0.99 &  =  &  ',num2str(Q_090_err_out_W0),'  &  ',num2str(Q_095_err_out_W0),'  &  ',num2str(Q_099_err_out_W0),'  &']);
disp(['W err outsample quantiles:&0.01 &0.05 & 0.10 &  =  &  ',num2str(Q_001_err_out_W0),'  &  ',num2str(Q_005_err_out_W0),'  &  ',num2str(Q_010_err_out_W0),'  &']);

disp(['WW1 err outsample quantiles:&0.90 &0.95 & 0.99 &  =  &  ',num2str(Q_090_err_out_WW10),'  &  ',num2str(Q_095_err_out_WW10),'  &  ',num2str(Q_099_err_out_WW10),'  &']);
disp(['WW1 err outsample quantiles:&0.01 &0.05 & 0.10 &  =  &  ',num2str(Q_001_err_out_WW10),'  &  ',num2str(Q_005_err_out_WW10),'  &  ',num2str(Q_010_err_out_WW10),'  &']);

disp(['WW2 err outsample quantiles:&0.90 &0.95 & 0.99 &  =  &  ',num2str(Q_090_err_out_WW20),'  &  ',num2str(Q_095_err_out_WW20),'  &  ',num2str(Q_099_err_out_WW20),'  &']);
disp(['WW2 err outsample quantiles:&0.01 &0.05 & 0.10 &  =  &  ',num2str(Q_001_err_out_WW20),'  &  ',num2str(Q_005_err_out_WW20),'  &  ',num2str(Q_010_err_out_WW20),'  &']);

disp(['WW3 err outsample quantiles:&0.90 &0.95 & 0.99 &  =  &  ',num2str(Q_090_err_out_WW30),'  &  ',num2str(Q_095_err_out_WW30),'  &  ',num2str(Q_099_err_out_WW30),'  &']);
disp(['WW3 err outsample quantiles:&0.01 &0.05 & 0.10 &  =  &  ',num2str(Q_001_err_out_WW30),'  &  ',num2str(Q_005_err_out_WW30),'  &  ',num2str(Q_010_err_out_WW30),'  &']);

disp('End Equal simulation set')

disp([' mean of N_next true values to be forecasted: ',num2str(mean(NextVals))]);
disp(['fails_detected_MLE=',num2str(fails_detected_MLE),' i.e. ',num2str(100*fails_detected_MLE/M),'% of M'])
disp(['fails_detected_W=',num2str(fails_detected_W),' i.e. ',num2str(100*fails_detected_W/M),'% of M'])
disp(['fails_detected_WW1=',num2str(fails_detected_WW1),' i.e. ',num2str(100*fails_detected_WW1/M),'% of M'])
disp(['fails_detected_WW2=',num2str(fails_detected_WW2),' i.e. ',num2str(100*fails_detected_WW2/M),'% of M'])
disp(['fails_detected_WW3=',num2str(fails_detected_WW3),' i.e. ',num2str(100*fails_detected_WW3/M),'% of M'])

disp(['NaNs_detected_MLE=',num2str(N_NaN_MLE),' i.e. ',num2str(100*N_NaN_MLE/M),'% of M'])
disp(['NaNs_detected_W=',num2str(N_NaN_W),' i.e. ',num2str(100*N_NaN_W/M),'% of M'])
disp(['NaNs_detected_WW1=',num2str(N_NaN_WW1),' i.e. ',num2str(100*N_NaN_WW1/M),'% of M'])
disp(['NaNs_detected_WW2=',num2str(N_NaN_WW2),' i.e. ',num2str(100*N_NaN_WW2/M),'% of M'])
disp(['NaNs_detected_WW3=',num2str(N_NaN_WW3),' i.e. ',num2str(100*N_NaN_WW3/M),'% of M'])

disp(['Total valid=',num2str(N_valid_all),' i.e. ',num2str(100*N_valid_all/M),'% of M'])

disp(['scenario= ',int2str(scenario)]);
disp([' method1: ',method1,', method2: ',method2,', method3: ',method3])
disp([' WW1: ',wavelet1,', WW2: ',wavelet2,', WW3: ',wavelet3])
disp([' level1: ',num2str(level1),', level2: ',num2str(level2),', level3: ',num2str(level3)])

disp('Check forecasting results above')

disp(['JumpsSim: ',JumpsSim])
disp(['JumpsEst: ',JumpsEst])

diary off;

end
