function [ average_x_values_jump, std_x_values_jump, upper_bound_jump, lower_bound_jump, x_values_jump ] = average_boom_bust_quinn( num_var, num_jump, jump_matrix, x )


% PURPOSE: CALCULATE AVERAGE AND TWO STANDARD DEVIATION INTERVAL IN BOOM/BUST
% EPISODES IN THE PERIOD OF THE BOOM/BUST ONE AND TWO PERIODS BEFORE AND ONE 
% AND TWO PERIODS AFTER THE EPISODES WHEN WE USE QUINN INDEX
% 
% INPUT:
%     -num_var: NUMBER OF TIME SERIES FOR THE VARIABLE FOR WHICH WE CALCULATE THE AVERAGE
%     -num_jump: NUMBER OF BOOM OR BUST
%     -jump_matrix: MATRIX WITH ALL BOOM/BUST EPISODES, YEAR AND COUNTRY WHERE THEY HAPPEN
%     -x: VARIABLE FOR WHICH WE CALCULATE THE AVERAGE AND INTERVAL
%     
% OUTPUT:
%     -average_x_values_jump: AVERAGE FOR THE VARIABLE IN CASE OF A TYPOLOGY OF EPISODE
%     -std_x_values_jump: STANDARD DEVIATION FOR THE VARIABLE IN CASE OF A TYPOLOGY OF EPISODE
%     -upper_bound_jump: UPPER BOUND OF THE INTERVAL CONFIDENCE FOR THE VARIABLE IN CASE OF A TYPOLOGY OF EPISODE
%     -lower_bound_jump: LOWER BOUND OF THE INTERVAL CONFIDENCE FOR THE VARIABLE IN CASE OF A TYPOLOGY OF EPISODE
%     -x_values_jump: MATRIX WITH ALL VALUES OF THE VARIABLE IN ALL CASES OF A TYPOLOGY OF EPISODES DIVIDED BY COUNTRY
%
% This function is needed because the Quinn Index is available only until 2007.
% So the available years are not 17, but 13 and the criteria has to be imposed 
% on the 12 and 13 observations


    %Create matrix for jumps (either booms or bustes)
    for j=1:num_var
        for i=1:num_jump
            if jump_matrix(i,2)==1
                x_values_jump (i,:,j)=[jump_matrix(i,1); nan; nan; x(jump_matrix(i,1),jump_matrix(i,2),j); x(jump_matrix(i,1),jump_matrix(i,2)+1,j); x(jump_matrix(i,1),jump_matrix(i,2)+2,j)];
            elseif jump_matrix(i,2)==2
                x_values_jump (i,:,j)=[jump_matrix(i,1); nan; x(jump_matrix(i,1),jump_matrix(i,2)-1,j); x(jump_matrix(i,1),jump_matrix(i,2),j); x(jump_matrix(i,1),jump_matrix(i,2)+1,j); x(jump_matrix(i,1),jump_matrix(i,2)+2,j)];
            elseif jump_matrix(i,2)==12
                x_values_jump (i,:,j)=[jump_matrix(i,1); x(jump_matrix(i,1),jump_matrix(i,2)-2,j); x(jump_matrix(i,1),jump_matrix(i,2)-1,j); x(jump_matrix(i,1),jump_matrix(i,2),j); x(jump_matrix(i,1),jump_matrix(i,2)+1,j); nan];    
            elseif jump_matrix(i,2)==13
                x_values_jump (i,:,j)=[jump_matrix(i,1); x(jump_matrix(i,1),jump_matrix(i,2)-2,j); x(jump_matrix(i,1),jump_matrix(i,2)-1,j); x(jump_matrix(i,1),jump_matrix(i,2),j); nan; nan];    
            else
                x_values_jump (i,:,j)=[jump_matrix(i,1); x(jump_matrix(i,1),jump_matrix(i,2)-2,j); x(jump_matrix(i,1),jump_matrix(i,2)-1,j); x(jump_matrix(i,1),jump_matrix(i,2),j); x(jump_matrix(i,1),jump_matrix(i,2)+1,j); x(jump_matrix(i,1),jump_matrix(i,2)+2,j)];
            end
        end
    end
 
    %Average for the different time horizons
    for j=1:num_var
        for i=1:5
            average_x_values_jump(1,i,j)=nanmean(x_values_jump(:,i+1,j),1);
        end
    end    
 
    %Standard Deviation
    for j=1:num_var
        for i=1:5
            std_x_values_jump(1,i,j)=nanstd(x_values_jump(:,i+1,j),0,1);
        end
    end
 
    %Confidence interval
    for j=1:num_var
        for i=1:5
            upper_bound_jump(1,i,j)=average_x_values_jump(1,i,j)+ 2*std_x_values_jump(1,i,j);
            lower_bound_jump(1,i,j)=average_x_values_jump(1,i,j)- 2*std_x_values_jump(1,i,j);
        end
    end

end

