
%Paper "Heterogeneity and the Macroeconomic Effects of Changes in
%Loan-to-Value Limits" by Emmanuel De Veirman and Jasper de Jong
%Journal of Money, Credit and Banking

addpath(genpath('JMCB_replic_deveirman-jong/MatlabTools'))

%Generate Pearson family distributions which are proxies for the
%Dutch empirical pre-cap distribution. Next, vary the standard
%deviation within the Pearson system.
%See Section 5 of the paper.

clear; close all; clc;

tic

%Directory to get data from
seed = 517859392;
rng(seed); %Seed for Random number generator

%This code requires a cross-sectional loan-to-value distribution.
%Because the Loan Level Data used in the paper are confidential,
%the LTV distributions are omitted.

%Define pre-cap distribution
nedpre = original_ltv_1011M8M7;

cd('JMCB_replic_deveirman-jong/Matlab/Newmatlabprograms')

%Compute moments of pre-cap distribution, which I will need in any case
mupre=mean(nedpre);
stpre=std(nedpre);
skpre=skewness(nedpre);
kupre=kurtosis(nedpre);

nedpremom=[mupre stpre skpre kupre]'; %collect moments of empirical pre-cap distribution,
%from top to bottom: mean stdev skew kurtosis

%FIRST: generate proxy for actual distribution

%Generate moments of one non-truncated Pearson distribution for which the
%truncated version has moments that are as close as possible to the
%moments of the empirical truncated distribution.
%For both the Pearson and the empirical distribution, "truncated" means:
%truncated to LTVs in (2,200].

%For each of first 4 moments: compute the difference between the actual and 
%a try-out value. %Minimize the absolute sum of these differences

fac1 = 1; %Proxy for the actual empirical distribution

obs = 1e6; %Number of observations in truncated Pearson family distribution

distyp = 1; %Distribution type. 1: Pearson. 2: normal.

if distyp == 1 %Pearson
    optmom = nan(4,size(fac1,2)); %Optimized moments
elseif distyp == 2 %Normal: only use population mean and stdev to optimize.
    optmom = nan(2,size(fac1,2)); %Optimized moments.
end

fvalvec = nan(1,size(fac1,2)); %Optimized objective
k=1;
for i = fac1
    target = [mupre i*stpre skpre kupre]';
    if distyp == 1
        x0=[target]; %Use target values (for trunc) as initial values for the moments of non-trunc dist
    elseif distyp == 2
        x0 = [target(1:2)]; %With normal distribution, only search across mean and stdev of raw normal
        %(but match all first four moments of capped distribution)
    end
    obsfac=1.1; %Factor to produce more observations for non-truncated P dist, out of which we keep obs
    options = optimset('TolFun',1,'TolX',1e-2,'MaxIter',1e4*length(x0),'MaxFunEvals',1e4*length(x0));
    %iteration criterion %'TolFun',1,'TolX',1e-2,
    vhuge=1e15; %For setting sum to huge number when undesirable outcome (when violate skew-kurt condition)
    [optmom(:,k), fvalvec(k),exitflag,output] = fminsearch(@pearstrunc,x0,options,target,obs,obsfac,vhuge,distyp)
    %MINIMIZE SUM OF ABSOLUTE DEVIATIONS
    k=k+1;
end
toc

if size(optmom,2) == 1
    x=optmom; %Optimized moments vector for non-truncated dist corresponding to trunc dist targeting empirical moments
end

%NEXT: generate alternative standard deviations for non-truncated Pearson family dist.
%In the end, use untruncated ánd truncated distributions
obs=1e6;
bigobs=round(obsfac*obs);
fac2=[0.80 0.85 0.90 0.95 1];  %Factor so as to vary (shrink) the standard deviation
Ptrunc=nan(obs,size(fac2,2)); %Matrix for truncated distributions from Pearson system
Puntrunc=nan(bigobs,size(fac2,2)); %Matrix with untruncated distributions from Pearson system
Puntrunc2=nan(obs,size(fac2,2)); %Matrix with untruncated distributions, selecting obs observations
Ptruncmom=nan(4,size(fac2,2)); %Matrix with moments for these distributions
Puntruncmom=nan(4,size(fac2,2)); %Matrix with moments for underlying non-truncated distributions
Puntrunc2mom=nan(4,size(fac2,2)); %Matrix with moments for non-truncated dist., selecting obs observations
Puntrunctyp=nan(1,size(fac2,2)); %Vector with Pearson types for non-truncated distributions

k=1;
for i=fac2
    %Generate from Pearson family distribution
    if distyp == 1 %actual Pearson
        [Puntrunc(:,k),Puntrunctyp(k)]=pearsrnd(x(1),i*x(2),x(3),x(4),bigobs,1);
    elseif distyp == 2 %normal
        [Puntrunc(:,k),Puntrunctyp(k)]=pearsrnd(x(1),i*x(2),0,3,bigobs,1);
    end
    Puntruncmom(:,k)=[mean(Puntrunc(:,k)) std(Puntrunc(:,k)) skewness(Puntrunc(:,k)) kurtosis(Puntrunc(:,k))]';
    %Moments of the Untruncated Pearson

    %Truncate
    pearstr=nan(size(Puntrunc(:,k)));
    j=0;
    for i=1:size(Puntrunc(:,k),1)
        if Puntrunc(i,k)>2 & Puntrunc(i,k)<=200 %select observations we wish to keep.
            j=j+1;  %indicates the number of observations we keep
            pearstr(j)=Puntrunc(i,k); %Store only values above 2 and not exceeding 200
        end
    end
    if j<obs %number of truncated observations is smaller than we need
        error('Number of truncated observations is too low. Increase obsfac.')
    elseif j>=obs %there are enough (and in general, too many) observations
        Ptrunc(:,k)=pearstr(1:obs); %select first obs observations and assign to Ptrunc (OUTPUT)
    end
    Ptruncmom(:,k)=[mean(Ptrunc(:,k)) std(Ptrunc(:,k)) skewness(Ptrunc(:,k)) kurtosis(Ptrunc(:,k))]'; %Moments of the truncated Pearson

    k=k+1;
end

Puntrunc2=Puntrunc(1:obs,:); %Only obs observations, not obsfac*obs (OUTPUT)
for k=1:size(Puntrunc2,2)
    Puntrunc2mom(:,k)=[mean(Puntrunc2(:,k)) std(Puntrunc2(:,k)) skewness(Puntrunc2(:,k)) kurtosis(Puntrunc2(:,k))]';
end
  
cd('JMCB_replic_deveirman-jong/LTV distributions')
 
writematrix(Puntrunc2, 'Puntrunc2.csv'); %1e6 observations from untruncated distributions
writematrix(Ptrunc, 'Ptrunc.csv'); %1e6 observations from truncated distributions

 


