function [pdf_extrap_matrix] = time_distrib_historic_extrap(pdf,timespan_pdf,pdf_scale_factor,timesteps_pdf_scale_factor,timespan_out)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%
%%% Applying linear steretch to a pdf using specified historic scaling
%%% factors 
%%%
%%% Output: extrapolated (stretched/scaled) pdf in every required
%%% historic output year (timespan_out) recorded as a matrix
%%%
%%% Dmitry Yumashev, 22/12/2019 
%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% make sure the time dimensions are correct to allow for matrix
%%% operations
if size(timespan_out,1) > 1 && size(timespan_out,2) == 1 % if the vector variable is a column
    timespan_out = transpose(timespan_out); % make it a row
end
%
if size(pdf,1) == 1 && size(pdf,2) > 1 % if the vector variable is a row
    pdf = transpose(pdf); % make it a column
end


%%% general settings

total_timesteps_out = numel(timespan_out);

total_timesteps_pdf = numel(timespan_pdf);

total_timesteps_scale_factor = numel(timesteps_pdf_scale_factor);

% pre-allocate memory for the output matrix; make sure it's set to zero to
% begin with!

pdf_extrap_matrix = zeros(total_timesteps_out,total_timesteps_out); % 1st argument: pdf time variable; 2nd argument: output year

%%% the scale factor is defined over longer time intervals, so that for
%%% time from timespan_scale_factor(j) to timespan_scale_factor(j+1) its
%%% value is scale_factor(j); assign scale factor values accordingly to
%%% every output point

scale_factor_extrap = ones(1,total_timesteps_out); % default value = 1 if the scale factor timespan does not cover all the output years, particularly in the near-term future forecast

for i=1:total_timesteps_out

    for j = 1:total_timesteps_scale_factor-1 % NOTE the last value       
    
        if timespan_out(i) >= timesteps_pdf_scale_factor(j) && timespan_out(i) < timesteps_pdf_scale_factor(j+1)

            scale_factor_extrap(i) = pdf_scale_factor(j);
            
            break % exit the j loop and move to the next step in the i loop

        end
        
    end
    
end


%%% now reconstruct cdf of the input pdf -- it will be used to carry out the
%%% transformation

cdf = zeros(total_timesteps_pdf,1); % use the same number of steps as pdf (not +1 ) for simplicity; make it a column

cdf(1) = 0; % make sure cdf starts from zero

for i=2:total_timesteps_pdf % NOTE the first value

    cdf(i) = cdf(i-1) + pdf(i-1) * (timespan_pdf(i) - timespan_pdf(i-1)); % use left-hand-side approximation (dafault method throughout)
    
end

% no need to worry about cdf(end) being set to 1 at this stage as it will
% be adjusted anyway after stretching


%%% now loop over output years and stretch the cdf according to the scale
%%% factor in each year

for j=1:total_timesteps_out

    new_pdf_time_limit = 1 + floor(scale_factor_extrap(j) * timespan_pdf(end));
    
    timespan_pdf_scaled = transpose(timespan_pdf(1):new_pdf_time_limit); % start with the first pdf timespan value and finish with the newly obtained last one; unit timestep 
    
    total_timesteps_pdf_scaled = numel(timespan_pdf_scaled);
    
    % allocate memory to the scaled cdf
    cdf_scaled = zeros(total_timesteps_pdf_scaled,1);
    pdf_scaled = zeros(total_timesteps_pdf_scaled,1); % same number of steps for simplicity
    
    % interpolate the values of the scaled cdf based on the surrounding steps in the original cdf
    
    % to do so, we need to map the new steps back onto the original timescale; note that the query steps won't always be integer
    timespan_pdf_original_query = timespan_pdf_scaled ./ scale_factor_extrap(j); % use scale factor for the current output step
    
    cdf_scaled = interp1(timespan_pdf, cdf, timespan_pdf_original_query, 'linear','extrap'); 
           
    %%% now we have got the new cdf and need to make sure it is normalised
    %%% to 1
    
    cdf_scaled = cdf_scaled / cdf_scaled(end);
    
    %%% now work out the scaled pdf 
    
    for i=1:total_timesteps_pdf_scaled-1 % NOTE the final step
    
        pdf_scaled(i) = (cdf_scaled(i+1) - cdf_scaled(i)) / (timespan_pdf_scaled(i+1) - timespan_pdf_scaled(i));
        
    end
    
    %%% record the scaled pdf for the given output year into output matrix;
    %%% cur the scaled pdf short if its timespan exceeds the output
    %%% timespan (no point having longer pdfs)
        
    if total_timesteps_pdf_scaled < total_timesteps_out
    
        pdf_extrap_matrix(1:total_timesteps_pdf_scaled,j) = pdf_scaled(:);
        pdf_extrap_matrix((total_timesteps_pdf_scaled+1):total_timesteps_out,j)= 0; % zeros in the remaining entries
    
    else % cut the pdf tail that falls beyond the output timespan
        
        pdf_extrap_matrix(:,j)= pdf_scaled(1:total_timesteps_out);
        
    end
    
end % close the loop over output years

%%%

end