function [mat_quant, matm, num_quant, ylimit]= fanChart(fore, hist,ylimit, datefinal,datestep)

% This function generates a fan chart. A fan chart aims to visualize the
% uncertainty that surrounds a sequence of point forecasts.
%
% Author: Marco Buchmann (mb-econ.net), ECB / Frankfurt University
% Date: January 2010
%
% Required input:
% -- a HxD matrix 'fore': the density forecasts D along the horizon 1:H
% -- a Tx1 vector 'hist': a vector of historical data
%
% Optional input:
% -- a string 'datefinal': referring to the end of the forecast horizon
% -- a scalar 'datestep': if e.g. set to 3, then labels quarterly 
%
% Output:
% -- plot

% Check input
if nargin<2
    error('Not enough input arguments')
end

% Parameters
H=size(fore,1);
T=size(hist,1);
Ttot=T+H;

%datefinal='Jun2009'; % date for which last forecast is available

ftitle='Out-of-sample forecast'; % title for plot
targetvar='Y'; % label for horizontal axis

% Date vector
if exist('datefinal')==1 %#ok<EXIST>
    if exist('datestep')==0 %#ok<EXIST>
       datestep=1; 
    end
    dvec=createdatevec(datefinal,Ttot,'backward');
    tmp_mnump=-(Ttot-1);
    % Generate date vector (displacement b/w start and end of test period in months)
    tmp_b=datevec((addtodate((addtodate(datenum(dvec(end,1),'mmmyyyy'),tmp_mnump,'month')),1,'day')));
    if tmp_b(1,2)==1 || tmp_b(1,2)==4 || tmp_b(1,2)==7 || tmp_b(1,2)==10
        sg=3:datestep:size(dvec,1);
    elseif tmp_b(1,2)==2 || tmp_b(1,2)==5 || tmp_b(1,2)==8 || tmp_b(1,2)==11
        sg=2:datestep:size(dvec,1);
    else
        sg=1:datestep:size(dvec,1);
    end
end

% Set relevant quantiles (in percent)
quantiles=[5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95];
%quantiles=[10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90];
%quantiles=[20 25 30 35 40 45 50 55 60 65 70 75 80] ;
num_quant=size(quantiles,2);

% Compute quantiles
mat_quant=NaN(H,num_quant);

for h=1:H % loop over horizon
    for q=1:num_quant % loop over quantiles
        mat_quant(h,q)=quantile(fore(h,:),quantiles(1,q)/100);
    end
end

% Add hist. observations to the plot matrix
 addv=[10^(-42) 10^(-43) 10^(-44) 10^(-45) 10^(-46) 10^(-47) 10^(-48) 10^(-49) 10^(-50) 0 10^(-50) 10^(-49) 10^(-48) 10^(-47) 10^(-46) 10^(-45) 10^(-44) 10^(-43) 10^(-42)];
%  addv=[10^(-43) 10^(-44) 10^(-45) 10^(-46) 10^(-47) 10^(-48) 10^(-49) 10^(-50) 0 10^(-50) 10^(-49) 10^(-48) 10^(-47) 10^(-46) 10^(-45) 10^(-44) 10^(-43)];
% addv=[10^(-45) 10^(-46) 10^(-47) 10^(-48) 10^(-49) 10^(-50) 0 10^(-50) 10^(-49) 10^(-48) 10^(-47) 10^(-46) 10^(-45) ];
% hm=hist;
% for i=1:num_quant-1
%     hm=[hm hist]; %#ok<AGROW>
% end
% mat_quant=[hm;mat_quant];

for t=1:T
   mat_quant(t,:)=mat_quant(t,:)+addv; 
end

% Prepare plot matrix for its use with the area function
matm=mat_quant;

for i=2:size(matm,2)
    matm(:,i)=matm(:,i)-mat_quant(:,i-1);
end
