%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Algorithm to estimate loan origination based on Adalid and Falagiarda (2019)
%   "How repayments manipulate our perceptions about loan dynamics after a boom"
%   Journal of Economics and Statistics, forthcoming
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all
close all
clc

% Dataset
[DATA,TEXT] = xlsread('Data_EuroArea','data1980');
Time = (DATA(:,1)+datenum('30/dec/1899'));
Dates = datevec(Time);
Dates(:,4:6) = [];

% Define variables, create matrices and set parameters
S = DATA(:,2); % Stock
mgr = DATA(:,3); % Monthly growth rate of the stock
N = DATA(:,4); % New business volumes
r = DATA(:,5); % Nominal interest rate
tm = DATA(:,6); % Maturity in months

Si = S(3,1); % Initial stock
ra = r./100; % Annual interest rate
rm = ra./12; % Monthly interest rate
T = size(N,1); % number of periods/loans
RP = zeros(tm(1,1)+1,1);
R = zeros(tm(1,1),1);
A = zeros(T+380+1,T+1); % Create matrix of repayments
rowsA = size(A,1);
colsA = size(A,2);

s = 0:0.01:1;
diff = zeros(length(s),1);
N_new = zeros(size(N,1),length(s));
Nn = zeros(rowsA,length(s));
S = cat(1,Si,zeros(rowsA,1));
grS = zeros(rowsA-11,1);
temp = zeros(12,1);
alpha = zeros(T,1);
w=4;

tic

% Start algorithm
progressbar
for i = 4:T
    for k = 1:length(s)
        
        N_new(w:end,k) = s(k).*N(w:end);
        RP(1,1) = N_new(i,k);
        
        for j = 1:tm(i,1)
            RP(j+1,:) = RP(1,1)*(1-(1+rm(i,1))^(j-tm(i,1)))/(1-(1+rm(i,1))^(-tm(i,1))); % Remaining principal (French loan formula)
            R(j,:) = RP(j,:)-RP(j+1,:); % Repayments
        end
        
        my_field = strcat('rp',num2str(i));
        rp.(my_field) =  RP;
        my_field2 = strcat('rep',num2str(i));
        rep.(my_field2) =  R;
        A(:,i+1) = cat(1,zeros(i,1),eval(['rep.rep' num2str(i)]),zeros(rowsA-size(R,1)-i,1));
        RP = zeros(tm(1,1)+1,1);
        R = zeros(tm(1,1),1);
        
        sumA = sum(A,2);
        date = datetime(Dates);
        
        for h = 1:colsA-1
            Nn(h,k) = N_new(h,k);
        end
        
        for x = 1:rowsA
            S(x+1) = S(x)+Nn(x,k)-sumA(x);
        end
        
        for z = 1:rowsA-11
            grS(z) = (S(z+1)/S(z)-1)*100; % monthly growth rate of stocks
        end
        
        for z = 1:rowsA-11
            agrS(z) = (S(z+12)/S(z)-1)*100; % annual growth rate of stocks
        end
        
        diff(k,1) = abs(mgr(i)-grS(i));
        if k == length(s)
            [M,I] = min(diff);
            N_new(w,:) = s(I).*N(w);
            alpha(i) = s(I);
            
            RP(1,1) = N_new(i,I);
            for j = 1:tm(i,1)
                RP(j+1,:) = RP(1,1)*(1-(1+rm(i,1))^(j-tm(i,1)))/(1-(1+rm(i,1))^(-tm(i,1))); % Remaining principal
                R(j,:) = RP(j,:)-RP(j+1,:); % Repayments
            end
            my_field = strcat('rp',num2str(i));
            rp.(my_field) =  RP;
            my_field2 = strcat('rep',num2str(i));
            rep.(my_field2) =  R;
            A(:,i+1) = cat(1,zeros(i,1),eval(['rep.rep' num2str(i)]),zeros(rowsA-size(R,1)-i,1));
            
            sumA = sum(A,2);
            date = datetime(Dates);
            
            Nn(i) = N_new(I);
            
            S(i+1) = S(i)+Nn(i)-sumA(i);
            
            grS(i) = (S(i+1)/S(i)-1)*100; % monthly growth rate of stocks
            
            agrS(i) = (S(i+12)/S(i)-1)*100; % annual growth rate of stocks
            
            w=w+1;
        end
        
    end
    progressbar(i/T)
end

toc

% Compute centered 7-month moving average of alpha
alpha_ma = movmean(alpha,7);

% Compute loan origination, monthly flows
LoanOrigination = N.*alpha_ma;

% Compute loan origination, annual flows
LoanOrigination_12m = movsum(LoanOrigination,[11 0]);

% Generate graph
figure('Color','white');
plot(date,LoanOrigination_12m,'-b');
legend('Loan origination, 12-month flows','Location','northwest');

