%This Matlab file is an example of the files used in 
%Maynard, A and Qiu, J., Public Insurance and Private Savings: 
%who is Affected and by How Much," Journal of Applied Econometrics, 
%forthcoming 

% This code calls several Matlab procedures written by Victor 
% Chernozhukov and Christian Hansen that can be dowloaded at:
% http://faculty.chicagogsb.edu/christian.hansen/research/iqrmat.zip
% This file itself is also based on a modification of the sample.m file
% included in the above mentioned zip file.

clear

% Load data

y = dlmread('../matlabdata/y.csv',',');
d = dlmread('../matlabdata/d.csv',',');
x = dlmread('../matlabdata/x.csv',',');
z = dlmread('../matlabdata/z.csv',',');

% y - dependent variable
% d - RHS endogenous variable
% x - covariates
% z - excluded instruments

w = [z,x,ones(size(y))];

dhat = w*inv(w'*w)*(w'*d);  % Generates linear projection of D on {X,Z} to use as instrument.
%disp('dhat');
%disp(dhat);
disp('Start Quantile Loop')
warning off;
%This makes the results more readible. Footnote 13
%of our paper provides further discussion.

for i = 0.1:0.1:0.9 %quantile index (i), start:increment:end
% Note: in some cases convergence was not obtained exactly at the
% specified decile. In this case results were re-estimated at
% quantiles 0.01 above or below the decile as shown for the median
% on the line below.
%for i = 0.49:0.02:0.51 %quantile index (i), start:increment:end
  disp('quantile i ');
  disp(i);
 
  disp('Next quantile regression');
  bqr50 = rq([d,x,ones(size(y))],y,i);          % Obtains QR estimates
  disp('Next quantile regression standard errors');
  bqr50 = vcqr(bqr50,y,[d,x],i)   % Obtains std errors and displays with estimates
            % First column - coefficients, Second column - Std errors
            % Estimates are displayed in order entered (i.e. d first, then x)
            % Intercept is always last.

   %A = (-500:10:1500)';        % Grid search region
  A=[]; %choose automatically
  disp('Next instrumental variables quantile regression estimates');            
  biqr50 = inv_qr(y,d,x,dhat,i,A);               % Obtains IQR estimates
  disp('Next instrumental variables quantile regression standard errors');            
  biqr50 = vciqr(biqr50,y,d,x,z,i)   % Obtains std errors and displays with estimates
            % First column - coefficients, Second column - Std errors
            % x estimates are displayed in order entered.
            % Treatment effect (d) is always first; intercept is always last.
  
end;
exit
