%% HOW TO USE THE PROGRAM

% In the Section 'SELECTED MODELS' (Lines 207 - 860) choose only one model.
% Currently, an aggregate model (i.e. which does not distinguish between
% loan categories (business, consumer etc) but uses the aggregate loan
% amounts, is selected (line 806). Only one of the models should be
% uncommented at each run of the programm.

%The models are categorised in the following sub-categories: 
% 1) 'Nominal models' (LHS = log(loans_bus_nf_short) OR log(loans_bus_nf_long) OR log(loans_cons) OR log(loans_mort) ); 
% 2) 'Models for break (pre- and post-crisis) analysis' (LHS = AS ABOVE); There are dummies for crisis and non-crisis periods 
% 3) 'Models with real variables' (LHS = log(loans_bus_nf_short_real) OR log(loans_bus_nf_long_real) OR log(loans_cons_real) OR log(loans_mort_real));
% 4) 'Other models' - others, various robustness checks etc

% When a model is chosen (by uncommenting it), then in order to present the
% graphs correctly, select (by uncommenting appropriately) the the corresponding number of variables in the
% code from lines 1320 to 1339. Please note that depending on the model you
% have selected, you should allow for one more variable, reflecting the
% constant term. Therefore, given the model selected in line 806, the lines
% 1333 and 1338 are uncommented, corresponding to 6 variables for the
% supply equation and 3 variables for the demand equation (including the
%constant term).

% Please note that due to confidentiality reasons, the data on
% non-performing loans are not included in the Excel file with the input data 
%(because these data have not been published by the Bank of Greece) therefore the
% respective models cannot be run. The respective commands in this file
% which process the non-performing loans data and run the models have been consequently commented in this
% file, so that the procedure followed remains transparent to the reader.


%% DEFINITION OF PROBLEM'S PARAMETERS
%--------------------------------------------------------------------------
        

    % Number of time observations - Period from 2003M12 until 2011M03 
    T=88;

    % Auxiliary parameter
    Infin=1000; % Just a proxy for infinity in the truncated normal
    
%% READ DATA
%--------------------------------------------------------------------------

file_data=['Data_credit_disequilibrium.xls'];
    % Read data
    
        % Macro real
        
        indprod_nsa = xlsread(file_data,1,'B2:B89');    
        une = xlsread(file_data,1,'C2:C89');    
        inf_dtk_yoy = xlsread(file_data,1,'D2:D89');    
        inf_endtk_yoy = xlsread(file_data,1,'E2:E89');    
        dtk_nsa = xlsread(file_data,1,'F2:F89');    
        fcapf_real_Nans = xlsread(file_data,1,'g2:g89');    
        gdp_real_Nans = xlsread(file_data,1,'h2:h89');        
        dummy_no_access = xlsread(file_data,1,'i2:i89');        
        
        gdp_nom_Nans = xlsread(file_data,1,'J2:J89');        
        
        indprod_1_nsa = xlsread(file_data,1,'o2:o89');    
        indprod_2_nsa = xlsread(file_data,1,'p2:p89');    
        
        % Time variable
        time = [1:1:88]';
    
        % Interest rates            
        ir_loan_outst = xlsread(file_data,2,'B2:B89');    	
        ir_loan_outst_hous = xlsread(file_data,2,'C2:C89');    		
        ir_loan_outst_firm = xlsread(file_data,2,'D2:D89');    		
        ir_loan_new = xlsread(file_data,2,'E2:E89');    		
        ir_loan_new_hous = xlsread(file_data,2,'F2:F89');    		
        ir_loan_new_firm = xlsread(file_data,2,'G2:G89');    		
        ir_margin = xlsread(file_data,2,'H2:H89');    	

        ir_loan_outst_firm_short = xlsread(file_data,2,'I2:I89');    		
        ir_loan_outst_firm_middle = xlsread(file_data,2,'J2:J89');    		
        ir_loan_outst_firm_long = xlsread(file_data,2,'K2:K89');    		
        
        ir_loan_new_mort = xlsread(file_data,2,'L2:L89');    		
        ir_loan_outst_mort = xlsread(file_data,2,'M2:M89');    		
        
        ir_loan_new_mort_1 = xlsread(file_data,2,'N2:N89');    		
        ir_loan_new_hous_1 = xlsread(file_data,2,'O2:O89');    		

        ir_loan_new_mort_2 = xlsread(file_data,2,'P2:P89');    		
        ir_loan_new_hous_2 = xlsread(file_data,2,'Q2:Q89');    		

        
        ir_loan_new_firm_short = xlsread(file_data,2,'R2:R89');    			
        ir_loan_new_firm_long = xlsread(file_data,2,'S2:S89');    		

        
        % Monetary policy
        ir_depf = xlsread(file_data,3,'B2:B89');    		
        ir_mlf = xlsread(file_data,3,'C2:C89');    		
        req_reserv = xlsread(file_data,3,'D2:D89');    	
        adj_euribor = xlsread(file_data,3,'e2:e89');    	
        
        % Credit deposits
        cre = xlsread(file_data,4,'B2:B89');    		 	 
        dep = xlsread(file_data,4,'C2:C89');    		 
        dep_1 = xlsread(file_data,4,'I2:I89');
        
        capital_ratio = xlsread(file_data,4,'W2:W89');

        %loans_priv = xlsread(file_data,4,'d2:d89');    		 
        loans_mort = xlsread(file_data,4,'e2:e89');    		 	
        loans_cons = xlsread(file_data,4,'f2:f89');    		 

        cre_bus_nf = xlsread(file_data,4,'G2:G89');    		 
        
        liabilities_ecb = xlsread(file_data,4,'h2:h89');    		 

        loans_bus_nf_short = xlsread(file_data,4,'r2:r89');    		 
        loans_bus_nf_long = xlsread(file_data,4,'s2:s89');    		 
        loans_bus_nf = xlsread(file_data,4,'t2:t89');    		 

        
        
        
        % Stock market
        ase = xlsread(file_data,5,'B2:B89');    		 	 	
        ase_ret = xlsread(file_data,5,'C2:C89');    		 	 	
        ase_vol_ma = xlsread(file_data,5,'D2:D89');    		 	 

        % Bond market
        yield_ger = xlsread(file_data,6,'B2:B89');    		 	 		
        yield_gr = xlsread(file_data,6,'C2:C89');    		 	 		
        spread = xlsread(file_data,6,'D2:D89');    		 	 	
        
        % NPL
%         npl_cons_Nans = xlsread(file_data,7,'B2:B89');    		 	 			
%         npl_mort_Nans = xlsread(file_data,7,'C2:C89');    		 	 			
%         npl_bus_Nans = xlsread(file_data,7,'D2:D89');    		 	 			
%         npl_Nans = xlsread(file_data,7,'E2:E89');    		 	 			
%         npl_rat_cons_Nans = xlsread(file_data,7,'F2:F89');    		 	 			
%         npl_rat_mort_Nans = xlsread(file_data,7,'G2:G89');    		 	 			
%         npl_rat_bus_Nans = xlsread(file_data,7,'H2:H89');    		 	 			
%         npl_rat_Nans = xlsread(file_data,7,'I2:I89');    		 	 		

        time1=[1:3:88];
        time2=[1:88];
        
        
        ii=find(not(isnan(fcapf_real_Nans)));
        fcapf_real_Nans=fcapf_real_Nans(ii);
        
        ii=find(not(isnan(gdp_real_Nans)));
        gdp_real_Nans=gdp_real_Nans(ii);
        
        ii=find(not(isnan(gdp_nom_Nans)));
        gdp_nom_Nans=gdp_nom_Nans(ii);
        
        
        
        % Interpolate variables that contain Nans
            
            fcapf_real = interp1(time1,fcapf_real_Nans,time2)';
            gdp_real = interp1(time1,gdp_real_Nans,time2)';
            gdp_nom = interp1(time1,gdp_nom_Nans,time2)';
            

        % Define a dummy which is 1 in case of negative GDP growth
        dummy_GDP_neg=zeros(T,1);
        dummy_GDP_neg(1)=0;
        for ij=2:T
            if (gdp_real(ij)-gdp_real(ij-1))>0
                dummy_GDP_neg(ij)=0;
            else
                dummy_GDP_neg(ij)=1;
            end;
        end;
        

      % Real quantities
      loans_bus_nf_short_real=loans_bus_nf_short.*dtk_nsa./100;
      loans_bus_nf_long_real=loans_bus_nf_long.*dtk_nsa./100;
      loans_mort_real=loans_mort.*dtk_nsa./100;
      loans_cons_real=loans_cons.*dtk_nsa./100;
      
      
      ir_loan_new_mort_real=ir_loan_new_mort-inf_dtk_yoy;
      ir_loan_new_hous_real=ir_loan_new_hous-inf_dtk_yoy;
      ir_loan_new_firm_short_real=ir_loan_new_firm_short-inf_dtk_yoy;
      ir_loan_new_firm_long_real=ir_loan_new_firm_long-inf_dtk_yoy;
      
      
        
    % Dependent variable    
    cre1000=log(cre/1000);
        
        
    % Variables for supply equation
    lendcap=dep/1000-req_reserv/1000;
    lendcap_1=dep_1/1000-req_reserv/1000;
    lendcap_real=lendcap.*dtk_nsa./100;
    % ir_loan_new
    % time
    
    ecb_dependence=liabilities_ecb./(dep);
    
    
    % Variables for demand equation
    % ir_loan_new
    % time
    ase1000=ase/1000;
    

% --------------------------------------------------------------------------      
      
%              SELECTED MODELS


%%             $$ Nominal models
%             x_s=[log(ase) log(lendcap)]; 
%             x_d=[(ir_loan_outst_firm_short) log(gdp_nom)];
%             LHS = log(loans_bus_nf_short);
%             Ns=2;
%             Nd=2;
%             supply_restrictions=[0; 0; 0;]; % 
%             demand_restrictions=[0; 0; 0;];

%              $$ Business long-term Loans
%              x_s=[log(lendcap) (ir_loan_new_firm_long) (ir_loan_new_firm_long).^2 (npl_rat_bus./npl_rat_cons)];
%              x_d=[(ir_loan_new_firm_long) log(gdp_nom)];
%              LHS = log(loans_bus_nf_long); 
%              Ns=4;
%              Nd=2;
%              supply_restrictions=[0; 0; 0; 0; 0;]; 
%              demand_restrictions=[0; 0; 0;];

%              $$ Consumer Loans
%              x_s=[log(lendcap)];
%              x_d=[(ir_loan_new_hous) log(gdp_nom)];
%              LHS = log(loans_cons);    
%              Ns=1;
%              Nd=2;
%              supply_restrictions=[0; 0;]; % 
%              demand_restrictions=[0; 0; 0;];
     
%              $$ Mortgage Loans
%              x_s=[log(lendcap)];
%              x_d=[(ir_loan_new_mort) log(gdp_nom)];
%              LHS = log(loans_mort);      
%              Ns=1;
%              Nd=2;
%              supply_restrictions=[0; 0;]; % 
%              demand_restrictions=[0; 0; 0;];

%% Models for break (pre- and post-crisis) analysis

%             x_s=[dummy_GDP_neg log(ase) log(lendcap) log(lendcap).*dummy_GDP_neg]; 
%             x_d=[dummy_GDP_neg (ir_loan_outst_firm_short) log(gdp_nom) log(gdp_nom).*dummy_GDP_neg];
%             LHS = log(loans_bus_nf_short);
%             Ns=2+2;
%             Nd=2+2;
%             supply_restrictions=[0; 0; 0; 0; 0;]; % 
%             demand_restrictions=[0; 0; 0; 0; 0;];

%              $$ Business long-term Loans
%              x_s=[dummy_GDP_neg log(lendcap) log(lendcap).*dummy_GDP_neg (ir_loan_new_firm_long) (ir_loan_new_firm_long).^2 (npl_rat_bus./npl_rat_cons)];
%              x_d=[dummy_GDP_neg (ir_loan_new_firm_long) log(gdp_nom) log(gdp_nom).*dummy_GDP_neg];
%              LHS = log(loans_bus_nf_long); 
%              Ns=4+2;
%              Nd=2+2;
%              supply_restrictions=[0; 0; 0; 0; 0; 0; 0;]; 
%              demand_restrictions=[0; 0; 0; 0; 0;];

%              $$ Consumer Loans
%              x_s=[dummy_GDP_neg log(lendcap) log(lendcap).*dummy_GDP_neg];
%              x_d=[dummy_GDP_neg (ir_loan_new_hous) log(gdp_nom) log(gdp_nom).*dummy_GDP_neg];
%              LHS = log(loans_cons);    
%              Ns=1+2;
%              Nd=2+2;
%              supply_restrictions=[0; 0; 0; 0;]; % 
%              demand_restrictions=[0; 0; 0; 0; 0;];
     
%              $$ Mortgage Loans
%              x_s=[dummy_GDP_neg log(lendcap) log(lendcap).*dummy_GDP_neg];
%              x_d=[dummy_GDP_neg (ir_loan_new_mort) log(gdp_nom) log(gdp_nom).*dummy_GDP_neg];
%              LHS = log(loans_mort);      
%              Ns=1+2;
%              Nd=2+2;
%              supply_restrictions=[0; 0; 0; 0;]; % 
%              demand_restrictions=[0; 0; 0; 0; 0;];


%%           Models with real variables  
%               $$ Business short-term Loans
%             x_s=[log(ase) log(lendcap_real)]; 
%             x_d=[(ir_loan_new_firm_short_real) log(gdp_real)];
%             LHS = log(loans_bus_nf_short_real);
%             Ns=2;
%             Nd=2;
%             supply_restrictions=[0; 0; 0;]; % 
%             demand_restrictions=[0; 0; 0;];

%              $$ Business long-term Loans
%              x_s=[log(lendcap_real) (ir_loan_new_firm_long_real) (ir_loan_new_firm_long_real).^2 (npl_rat_bus_real./npl_rat_cons_real)];
%              x_d=[(ir_loan_new_firm_long_real) log(gdp_real)];
%              LHS = log(loans_bus_nf_long_real); 
%              Ns=4;
%              Nd=2;
%              supply_restrictions=[0; 0; 0; 0; 0;]; 
%              demand_restrictions=[0; 0; 0;];

%              $$ Consumer Loans
%              x_s=[log(lendcap_real)];
%              x_d=[(ir_loan_new_hous_real) log(gdp_real)];
%              LHS = log(loans_cons_real);    
%              Ns=1;
%              Nd=2;
%              supply_restrictions=[0; 0;]; % 
%              demand_restrictions=[0; 0; 0;];
     
%              $$ Mortgage Loans
%              x_s=[log(lendcap_real)];
%              x_d=[(ir_loan_new_mort) log(gdp_real)];
%              LHS = log(loans_mort_real);      
%              Ns=1;
%              Nd=2;
%              supply_restrictions=[0; 0;]; % 
%              demand_restrictions=[0; 0; 0;];


%% OTHER MODELS
% --------------------------------------------------------------------------      

%              MODEL 1_1       
%              x_s=[log(lendcap)];
%              x_d=[(ir_loan_new_hous) log(gdp_nom)];
%              LHS = log(loans_cons);    
%              Ns=1;
%              Nd=2;
%              supply_restrictions=[0; 0;]; 
%              demand_restrictions=[0; 0; 0;];

%              MODEL 1_1_lagged       
%              x_s=[log(lendcap_1)];
%              x_d=[(ir_loan_new_hous) log(gdp_nom)];
%              LHS = log(loans_cons);    
%              Ns=1;
%              Nd=2;
%              supply_restrictions=[0; 0;]; 
%              demand_restrictions=[0; 0; 0;];

     
%              MODEL 3_1
%              x_s=[log(lendcap)];
%              x_d=[(ir_loan_new_mort) log(gdp_nom)];
%              LHS = log(loans_mort);      
%              Ns=1;
%              Nd=2;
%              supply_restrictions=[0; 0;]; 
%              demand_restrictions=[0; 0; 0;];


%              MODEL 3_1_lagged
%              x_s=[log(lendcap_1)];
%              x_d=[(ir_loan_new_mort) log(gdp_nom)];
%              LHS = log(loans_mort);      
%              Ns=1;
%              Nd=2;
%              supply_restrictions=[0; 0;]; 
%              demand_restrictions=[0; 0; 0;];



%              MODEL 1_2      
%              x_s=[log(ase) log(lendcap)];
%              x_d=[(ir_loan_new_hous) log(gdp_nom)];
%              LHS = log(loans_cons);    
%              Ns=2;
%              Nd=2;
%              supply_restrictions=[0; 0; 0;]; 
%              demand_restrictions=[0; 0; 0;];

     
%              MODEL 3_2
%              x_s=[log(ase) log(lendcap)];
%              x_d=[(ir_loan_new_mort) log(gdp_nom)];
%              LHS = log(loans_mort);      
%              Ns=2;
%              Nd=2;
%              supply_restrictions=[0; 0; 0;]; 
%              demand_restrictions=[0; 0; 0;];



% --------------------------------------------------------------------------      
      
%              MODEL 1_3      
%              x_s=[log(lendcap) (ir_loan_new_hous) (ir_loan_new_hous).^2];
%              x_d=[(ir_loan_new_hous) log(gdp_nom)];
%              LHS = log(loans_cons);    
%              Ns=3;
%              Nd=2;
%              supply_restrictions=[0; 0; 0; 0;]; 
%              demand_restrictions=[0; 0; 0;];

     
%              MODEL 3_3
%              x_s=[log(lendcap) (ir_loan_new_mort) (ir_loan_new_mort).^2];
%              x_d=[(ir_loan_new_mort) log(gdp_nom)];
%              LHS = log(loans_mort);      
%              Ns=3;
%              Nd=2;
%              supply_restrictions=[0; 0; 0; 0;]; 
%              demand_restrictions=[0; 0; 0;];

      
% --------------------------------------------------------------------------      
      
%              MODEL 1_4      
%              x_s=[log(lendcap) (ir_loan_new_hous) (ir_loan_new_hous).^2 (npl_rat_cons) (npl_rat_bus)];
%              x_d=[(ir_loan_new_hous) log(gdp_nom)];
%              LHS = log(loans_cons);    
%              Ns=5;
%              Nd=2;
%              supply_restrictions=[0; 0; 0; 0; 0; 0;]; 
%              demand_restrictions=[0; 0; 0;];

     
   % MODEL 3_4
%              x_s=[log(lendcap) (ir_loan_new_mort) (ir_loan_new_mort).^2 (npl_rat_mort) npl_rat_bus];
%              x_d=[(ir_loan_new_mort) log(gdp_nom)];
%              LHS = log(loans_mort);      
%              Ns=5;
%              Nd=2;
%              supply_restrictions=[0; 0; 0; 0; 0; 0;]; 
%              demand_restrictions=[0; 0; 0;];
% 

%-------------------

%             MODEL 1_5      
%             x_s=[log(ase) log(lendcap) (ir_loan_new_hous) (ir_loan_new_hous).^2];
%             x_d=[(ir_loan_new_hous) log(gdp_nom)];
%             LHS = log(loans_cons);
%             Ns=4;
%             Nd=2;
%             supply_restrictions=[0; 0; 0; 0; 0;]; 
%             demand_restrictions=[0; 0; 0;];

                % MODEL 3_5
%             x_s=[log(ase) log(lendcap) (ir_loan_new_mort) (ir_loan_new_mort).^2];
%             x_d=[(ir_loan_new_mort) log(gdp_nom)];
%             LHS = log(loans_mort);       
%             Ns=4;
%             Nd=2;
%             supply_restrictions=[0; 0; 0; 0; 0;]; 
%             demand_restrictions=[0; 0; 0;];

%-------------------




%------------------- 
                 % Model 2_1_short
%             x_s=[log(lendcap)]; 
%             x_d=[(ir_loan_outst_firm_short) log(gdp_nom)];
%             LHS = log(loans_bus_nf_short);
%             Ns=1;
%             Nd=2;
%             supply_restrictions=[0; 0;]; 
%             demand_restrictions=[0; 0; 0;];
    
                % Model 2_1_long    
%             x_s=[log(lendcap)]; 
%             x_d=[(ir_loan_outst_firm_long) log(gdp_nom)];
%             LHS = log(loans_bus_nf_long);
%             Ns=1;
%             Nd=2;
%             supply_restrictions=[0; 0;]; 
%             demand_restrictions=[0; 0; 0;];


%-------------------
%                 % Model 2_2_short
%             x_s=[log(ase) log(lendcap)]; 
%             x_d=[(ir_loan_outst_firm_short) log(gdp_nom)];
%             LHS = log(loans_bus_nf_short);
%             Ns=2;
%             Nd=2;
%             supply_restrictions=[0; 0; 0;]; 
%             demand_restrictions=[0; 0; 0;];
    
            % Model 2_2_short_lagged
%             x_s=[log(ase) log(lendcap_1)]; 
%             x_d=[(ir_loan_outst_firm_short) log(gdp_nom)];
%             LHS = log(loans_bus_nf_short);
%             Ns=2;
%             Nd=2;
%             supply_restrictions=[0; 0; 0;]; 
%             demand_restrictions=[0; 0; 0;];
    
            
            
                % Model 2_2_long    
%             x_s=[log(ase) log(lendcap)]; 
%             x_d=[(ir_loan_outst_firm_long) log(gdp_nom)];
%             LHS = log(loans_bus_nf_long);
%             Ns=2;
%             Nd=2;
%             supply_restrictions=[0; 0; 0;]; 
%             demand_restrictions=[0; 0; 0;];


%-------------------
                 % Model (2_3_short) - Short-term business with Interest rates 
%              x_s=[log(lendcap) (ir_loan_new_firm_short) (ir_loan_new_firm_short).^2];
%              x_d=[(ir_loan_new_firm_short) log(gdp_nom)];
%              LHS = log(loans_bus_nf_short);  
%              Ns=3;
%              Nd=2;
%              supply_restrictions=[0; 0; 0; 0;]; 
%              demand_restrictions=[0; 0; 0;];


                % Model (2_3_long) - Long-term business with Interest rates 
%             x_s=[log(lendcap) (ir_loan_new_firm_long) (ir_loan_new_firm_long).^2 ir_loan_new_firm_short];
%             x_d=[(ir_loan_new_firm_long) log(gdp_nom)];
%             LHS = log(loans_bus_nf_long);  
%              Ns=4;
%              Nd=2;
%              supply_restrictions=[0; 0; 0; 0; 0;]; 
%              demand_restrictions=[0; 0; 0;];


%-------------------


                 % MODEL  2_4_short
%              x_s=[log(lendcap) (ir_loan_new_firm_short) (ir_loan_new_firm_short).^2 (npl_rat_bus)];
%              x_d=[(ir_loan_new_firm_short) log(gdp_nom)];
%              LHS = log(loans_bus_nf_short);
%              Ns=4;
%              Nd=2;
%              supply_restrictions=[0; 0; 0; 0; 0;]; 
%              demand_restrictions=[0; 0; 0;];
     
                 %MODEL  2_4_long  
%              x_s=[log(lendcap) (ir_loan_new_firm_long) (ir_loan_new_firm_long).^2 (npl_rat_bus./npl_rat_cons)];
%              x_d=[(ir_loan_new_firm_long) log(gdp_nom)];
%              LHS = log(loans_bus_nf_long); 
%              Ns=4;
%              Nd=2;
%              supply_restrictions=[0; 0; 0; 0; 0;]; 
%              demand_restrictions=[0; 0; 0;];
   
                 %MODEL  2_4_long-lagged  
%              x_s=[log(lendcap_1) (ir_loan_new_firm_long) (ir_loan_new_firm_long).^2 (npl_rat_bus./npl_rat_cons)];
%              x_d=[(ir_loan_new_firm_long) log(gdp_nom)];
%              LHS = log(loans_bus_nf_long); 
%              Ns=4;
%              Nd=2;
%              supply_restrictions=[0; 0; 0; 0; 0;]; 
%              demand_restrictions=[0; 0; 0;];

%-------------------

                 % MODEL  2_5_short 
%              x_s=[log(ase) log(lendcap) (ir_loan_new_firm_short) (ir_loan_new_firm_short).^2];
%              x_d=[(ir_loan_new_firm_short) log(gdp_nom)];
%              LHS = log(loans_bus_nf_short); 
%              Ns=4;
%              Nd=2;
%              supply_restrictions=[0; 0; 0; 0; 0;]; 
%              demand_restrictions=[0; 0; 0;];

         
                 % MODEL  2_5_long  
%              x_s=[log(ase) log(lendcap) (ir_loan_new_firm_long) (ir_loan_new_firm_long).^2];
%              x_d=[(ir_loan_new_firm_long) log(gdp_nom)];
%              LHS = log(loans_bus_nf_long);
%              Ns=4;
%              Nd=2;
%              supply_restrictions=[0; 0; 0; 0; 0;]; 
%              demand_restrictions=[0; 0; 0;];

%-------------------

% ROBUSTNESS CHECK MODELS - PREFERRED MODELS with extra variables

%--------------------------------------------------------------------
% INCLUDING ecb_dependence

%              MODEL 1_6      
%              x_s=[log(lendcap) ecb_dependence];
%              x_d=[(ir_loan_new_hous) log(gdp_nom)];
%              LHS = log(loans_cons);    
%              Ns=2;
%              Nd=2;
%              supply_restrictions=[0; 0; 0;]; 
%              demand_restrictions=[0; 0; 0;];

%              MODEL 3_6
%              x_s=[log(lendcap) ecb_dependence];
%              x_d=[(ir_loan_new_mort) log(gdp_nom)];
%              LHS = log(loans_mort);      
%              Ns=2;
%              Nd=2;
%              supply_restrictions=[0; 0; 0;]; 
%              demand_restrictions=[0; 0; 0;];

                % Model 2_6_short
%             x_s=[log(ase) log(lendcap) ecb_dependence]; 
%             x_d=[(ir_loan_outst_firm_short) log(gdp_nom)];
%             LHS = log(loans_bus_nf_short);
%             Ns=3;
%             Nd=2;
%             supply_restrictions=[0; 0; 0; 0;]; 
%             demand_restrictions=[0; 0; 0;];


                 %MODEL  2_6_long  
%              x_s=[log(lendcap) (ir_loan_new_firm_long) (ir_loan_new_firm_long).^2 (npl_rat_bus./npl_rat_cons) ecb_dependence];
%              x_d=[(ir_loan_new_firm_long) log(gdp_nom)];
%              LHS = log(loans_bus_nf_long); 
%              Ns=5;
%              Nd=2;
%              supply_restrictions=[0; 0; 0; 0; 0; 0;]; 
%              demand_restrictions=[0; 0; 0;];

%---------------------

%--------------------------------------------------------------------
% INCLUDING ASE in demand

%              MODEL 1_7      
%              x_s=[log(lendcap)];
%              x_d=[(ir_loan_new_hous) log(gdp_nom) log(ase)];
%              LHS = log(loans_cons);    
%              Ns=1;
%              Nd=3;
%              supply_restrictions=[0; 0;]; 
%              demand_restrictions=[0; 0; 0; 0;];

%              MODEL 3_7
%              x_s=[log(lendcap)];
%              x_d=[(ir_loan_new_mort) log(gdp_nom) log(ase)];
%              LHS = log(loans_mort);      
%              Ns=1;
%              Nd=3;
%              supply_restrictions=[0; 0;]; 
%              demand_restrictions=[0; 0; 0; 0;];

                % Model 2_7_short
%             x_s=[log(ase) log(lendcap)]; 
%             x_d=[(ir_loan_outst_firm_short) log(gdp_nom) log(ase)];
%             LHS = log(loans_bus_nf_short);
%             Ns=2;
%             Nd=3;
%             supply_restrictions=[0; 0; 0;]; 
%             demand_restrictions=[0; 0; 0; 0;];


                 %MODEL  2_7_long  
%              x_s=[log(lendcap) (ir_loan_new_firm_long) (ir_loan_new_firm_long).^2 (npl_rat_bus./npl_rat_cons)];
%              x_d=[(ir_loan_new_firm_long) log(gdp_nom) log(ase)];
%              LHS = log(loans_bus_nf_long); 
%              Ns=4;
%              Nd=3;
%              supply_restrictions=[0; 0; 0; 0; 0;]; 
%              demand_restrictions=[0; 0; 0; 0;];

%---------------------

%  INCLUDING  inflation in both supply and demand



%              MODEL 1_8      
%              x_s=[log(lendcap) inf_dtk_yoy];
%              x_d=[(ir_loan_new_hous) log(gdp_nom) inf_dtk_yoy];
%              LHS = log(loans_cons);    
%              Ns=2;
%              Nd=3;
%              supply_restrictions=[0; 0; 0;]; 
%              demand_restrictions=[0; 0; 0; 0;];

%              MODEL 3_8
%              x_s=[log(lendcap) inf_dtk_yoy];
%              x_d=[(ir_loan_new_mort) log(gdp_nom) inf_dtk_yoy];
%              LHS = log(loans_mort);      
%              Ns=2;
%              Nd=3;
%              supply_restrictions=[0; 0; 0;]; 
%              demand_restrictions=[0; 0; 0; 0;];

                % Model 2_8_short
%             x_s=[log(ase) log(lendcap) inf_dtk_yoy]; 
%             x_d=[(ir_loan_outst_firm_short) log(gdp_nom) inf_dtk_yoy];
%             LHS = log(loans_bus_nf_short);
%             Ns=3;
%             Nd=3;
%             supply_restrictions=[0; 0; 0; 0;]; 
%             demand_restrictions=[0; 0; 0; 0;];


                 %MODEL  2_8_long  
%              x_s=[log(lendcap) (ir_loan_new_firm_long) (ir_loan_new_firm_long).^2 (npl_rat_bus./npl_rat_cons) inf_dtk_yoy];
%              x_d=[(ir_loan_new_firm_long) log(gdp_nom) log(ase) inf_dtk_yoy];
%              LHS = log(loans_bus_nf_long); 
%              Ns=5;
%              Nd=4;
%              supply_restrictions=[0; 0; 0; 0; 0; 0;]; 
%              demand_restrictions=[0; 0; 0; 0; 0;];

%-----------------------------------------------------
% AGGREGATE MODELS

loans_tot=loans_mort+loans_bus_nf_long+loans_cons+loans_bus_nf_short;
w_ir=(loans_mort./loans_tot).*ir_loan_new_mort + (loans_bus_nf_long./loans_tot).* ir_loan_new_firm_long + (loans_bus_nf_short./loans_tot).* ir_loan_new_firm_short + (loans_mort./loans_tot).* ir_loan_new_hous;
%w_npl=(loans_mort./loans_tot).*npl_rat_mort + ((loans_bus_nf_long+loans_bus_nf_short)./loans_tot).* npl_rat_bus + (loans_mort./loans_tot).* npl_rat_cons;

                 % Model Aggr_1
%             x_s=[log(lendcap)]; 
%             x_d=[(w_ir) log(gdp_nom)];
%             LHS = log(loans_tot);
%             Ns=1;
%             Nd=2;
%             supply_restrictions=[0; 0;]; 
%             demand_restrictions=[0; 0; 0;];

                % Model Aggr_2
%             x_s=[log(ase) log(lendcap)]; 
%             x_d=[w_ir log(gdp_nom)];
%             LHS = log(loans_tot);
%             Ns=2;
%             Nd=2;
%             supply_restrictions=[0; 0; 0;]; 
%             demand_restrictions=[0; 0; 0;];

                 % Model Aggr_3_short
%              x_s=[log(lendcap) w_ir w_ir.^2];
%              x_d=[w_ir log(gdp_nom)];
%              LHS = log(loans_tot);  
%              Ns=3;
%              Nd=2;
%              supply_restrictions=[0; 0; 0; 0;]; 
%              demand_restrictions=[0; 0; 0;];

                % Model Aggr_4_short
%              x_s=[log(lendcap) w_ir w_ir.^2 w_npl];
%              x_d=[w_ir log(gdp_nom)];
%              LHS = log(loans_tot);
%              Ns=4;
%              Nd=2;
%              supply_restrictions=[0; 0; 0; 0; 0;]; 
%              demand_restrictions=[0; 0; 0;];

                % Model Aggr_5_short
%              x_s=[log(ase) log(lendcap) w_ir w_ir.^2];
%              x_d=[w_ir log(gdp_nom)];
%              LHS = log(loans_tot); 
%              Ns=4;
%              Nd=2;
%              supply_restrictions=[0; 0; 0; 0; 0;]; 
%              demand_restrictions=[0; 0; 0;];

%----------------------------------




% AGGREGATE MODELS with capital_ratio

                 % Model Aggr_1
%             x_s=[log(lendcap) capital_ratio]; 
%             x_d=[(w_ir) log(gdp_nom)];
%             LHS = log(loans_tot);
%             Ns=2;
%             Nd=2;
%             supply_restrictions=[0; 0; 0;]; 
%             demand_restrictions=[0; 0; 0;];

                % Model Aggr_2
%             x_s=[log(ase) capital_ratio log(lendcap)]; 
%             x_d=[w_ir log(gdp_nom)];
%             LHS = log(loans_tot);
%             Ns=3;
%             Nd=2;
%             supply_restrictions=[0; 0; 0; 0;]; 
%             demand_restrictions=[0; 0; 0;];

                 % Model Aggr_3_short
%              x_s=[log(lendcap) capital_ratio w_ir w_ir.^2];
%              x_d=[w_ir log(gdp_nom)];
%              LHS = log(loans_tot);  
%              Ns=4;
%              Nd=2;
%              supply_restrictions=[0; 0; 0; 0; 0;]; 
%              demand_restrictions=[0; 0; 0;];

                % Model Aggr_4_short
%              x_s=[log(lendcap) capital_ratio w_ir w_ir.^2 w_npl];
%              x_d=[w_ir log(gdp_nom)];
%              LHS = log(loans_tot);
%              Ns=5;
%              Nd=2;
%              supply_restrictions=[0; 0; 0; 0; 0; 0;]; 
%              demand_restrictions=[0; 0; 0;];

                % Model Aggr_5_short
             x_s=[log(ase) capital_ratio log(lendcap) w_ir w_ir.^2];
             x_d=[w_ir log(gdp_nom)];
             LHS = log(loans_tot); 
             Ns=5;
             Nd=2;
             supply_restrictions=[0; 0; 0; 0; 0; 0;]; 
             demand_restrictions=[0; 0; 0;];


%-----------------------------------------------------
% AGGREGATE MODELS with yields

                 % Model Aggr_1
%             x_s=[log(lendcap)]; 
%             x_d=[(w_ir-yield_ger) log(gdp_nom)];
%             LHS = log(loans_tot);
%             Ns=1;
%             Nd=2;
%             supply_restrictions=[0; 0;]; 
%             demand_restrictions=[0; 0; 0;];

                % Model Aggr_2
%             x_s=[log(ase) log(lendcap)]; 
%             x_d=[w_ir-yield_ger log(gdp_nom)];
%             LHS = log(loans_tot);
%             Ns=2;
%             Nd=2;
%             supply_restrictions=[0; 0; 0;]; 
%             demand_restrictions=[0; 0; 0;];

                 % Model Aggr_3_short
%              x_s=[log(lendcap) w_ir w_ir.^2];
%              x_d=[w_ir-yield_ger log(gdp_nom)];
%              LHS = log(loans_tot);  
%              Ns=3;
%              Nd=2;
%              supply_restrictions=[0; 0; 0; 0;]; 
%              demand_restrictions=[0; 0; 0;];

                % Model Aggr_4_short
%              x_s=[log(lendcap) w_ir w_ir.^2 w_npl];
%              x_d=[w_ir-yield_ger log(gdp_nom)];
%              LHS = log(loans_tot);
%              Ns=4;
%              Nd=2;
%              supply_restrictions=[0; 0; 0; 0; 0;]; 
%              demand_restrictions=[0; 0; 0;];

                % Model Aggr_5_short
%              x_s=[log(ase) log(lendcap) w_ir w_ir.^2];
%              x_d=[w_ir-yield_ger log(gdp_nom)];
%              LHS = log(loans_tot); 
%              Ns=4;
%              Nd=2;
%              supply_restrictions=[0; 0; 0; 0; 0;]; 
%              demand_restrictions=[0; 0; 0;];

%----------------------------------

%% DEFINITION OF SOLUTION'S PARAMETERS
%--------------------------------------------------------------------------
        
    Coeffs=zeros(1000,Ns+1); %1000 is the number of iterations saved
    Coeffd=zeros(1000,Nd+1);
    
    stds_mat=zeros(1000,1); %1000 is the number of iterations saved
    stdd_mat=zeros(1000,1);

    ss2_mat=zeros(1000,1); %1000 is the number of iterations saved
    sd2_mat=zeros(1000,1);
    
    ys_mat=zeros(1000,T); %1000 is the number of iterations saved
    yd_mat=zeros(1000,T); 
    
    esupply_mat=zeros(1000,T); %1000 is the number of iterations saved
    edemand_mat=zeros(1000,T);
            
    lamda_mat=zeros(1000,T);
    ys_rel_yd=zeros(1000,T);
    
        %corresponding variables for the burnoff phase
        Coeffs_burnoff=zeros(1000,Ns+1); %1000 is the number of iterations saved
        Coeffd_burnoff=zeros(1000,Nd+1);
    
        stds_mat_burnoff=zeros(1000,1); %1000 is the number of iterations saved
        stdd_mat_burnoff=zeros(1000,1);

        ss2_mat_burnoff=zeros(1000,1); %1000 is the number of iterations saved
        sd2_mat_burnoff=zeros(1000,1);
    
        ys_mat_burnoff=zeros(1000,T); %1000 is the number of iterations saved
        yd_mat_burnoff=zeros(1000,T); 
    
        esupply_mat_burnoff=zeros(1000,T); %1000 is the number of iterations saved
        edemand_mat_burnoff=zeros(1000,T);
            
        lamda_mat_burnoff=zeros(1000,T);
    
    
    
    
%% OLS ESTIMATION (used for initialization and definition of posterior
%% parameters)
%--------------------------------------------------------------------------

    % Define LHS and LHS
    Constant_column=ones(T,1);
    
    ys=zeros(T,1);
    yd=zeros(T,1);
    
    RHS_s=[Constant_column x_s]; 
    RHS_d=[Constant_column x_d];

    [Coeffs_OLS_init,bint_s,r_s,rint_s,stats_s] = regress(LHS,RHS_s);
    [Coeffd_OLS_init,bint_d,r_d,rint_d,stats_d] = regress(LHS,RHS_d);

    ss2_OLS_init=((LHS-RHS_s*Coeffs_OLS_init)'*(LHS-RHS_s*Coeffs_OLS_init))/(T-Ns);
    sd2_OLS_init=((LHS-RHS_d*Coeffd_OLS_init)'*(LHS-RHS_d*Coeffd_OLS_init))/(T-Nd);
    
    figure;plot(RHS_d*Coeffd_OLS_init,'r');hold on;plot(RHS_s*Coeffs_OLS_init,'b');hold on;plot(LHS,'k');legend('PotDemand','PotSupply','Actual credit');
    title('OLS');
    
    Coeffs_OLS=Coeffs_OLS_init;
    Coeffd_OLS=Coeffd_OLS_init;
    

    
%% PRIOR HYPERPARAMETERS (I define prior parameters after the posterior because I use OLS estimations)

variance_factor=0.0001;

shift_factor=1;

% Coefficients
    Coeffs_prior=shift_factor*Coeffs_OLS_init;
    ii=find((sign(supply_restrictions)~=sign(Coeffs_OLS_init))&supply_restrictions~=0);
    Coeffs_prior(ii)=-shift_factor*Coeffs_prior(ii);
    
    
    Coeffd_prior=Coeffd_OLS_init;
    ii=find((sign(demand_restrictions)~=sign(Coeffd_OLS_init))&demand_restrictions~=0);
    Coeffd_prior(ii)=-Coeffd_prior(ii);
    
    
% Variance-covariance matrices   
    std_s_OLS_derived=sqrt(abs(0.5*bint_s(:,2)-Coeffs_OLS_init));
    Var_diagonal=variance_factor*std_s_OLS_derived.^2;
        ii=find((sign(supply_restrictions)~=sign(Coeffs_OLS_init))&supply_restrictions~=0);
        if not(isempty(ii))
            scaling_factor=variance_factor*std_s_OLS_derived.^2;
            Var_diagonal(ii)=Var_diagonal(ii).*scaling_factor;
        end;        
    Var_prior_s=diag(Var_diagonal);
    
    std_d_OLS_derived=sqrt(abs(0.5*bint_d(:,2)-Coeffd_OLS_init));
    Var_diagonal=variance_factor*std_d_OLS_derived.^2;
        ii=find((sign(demand_restrictions)~=sign(Coeffd_OLS_init))&demand_restrictions~=0);
        if not(isempty(ii))
            scaling_factor=variance_factor*std_d_OLS_derived.^2;
            Var_diagonal(ii)=Var_diagonal(ii).*scaling_factor;
        end;        
    Var_prior_d=diag(Var_diagonal);
    
% Degrees of freedom
    dofs_prior=T-Ns;
    dofd_prior=T-Nd;
    
% Standard deviations
    std_factor=1;
    ss2_prior=std_factor*stats_s(4);
    sd2_prior=std_factor*stats_d(4);
    
    
%% MCMC Algorithm
%--------------------------------------------------------------------------

    % Definition of variables used in MCMC algorithm
    lamda=zeros(T,1);

    % Some OLS parameters (except Coeffs_OLS and Coeffd_OLS)
    dofs_OLS=T-Ns;  %degrees of freedom for the supply equation
    dofd_OLS=T-Nd;  %degrees of freedom for the demand equation
    
    XXs=(RHS_s'*RHS_s);
    XXd=(RHS_d'*RHS_d);
    
    % Parameters of the posterior
    Vs_posterior=inv(inv(Var_prior_s)+RHS_s'*RHS_s);
    Vd_posterior=inv(inv(Var_prior_d)+RHS_d'*RHS_d);

    % Decompose the covariance matrix into the low triangular
    % Cholesky decomposition matrix
        VAR_s_mat=Vs_posterior;
    Chol_s_mat=chol(VAR_s_mat,'lower'); 

    VAR_d_mat=Vd_posterior;
    Chol_d_mat=chol(VAR_d_mat,'lower');

    %(Eq. 3.12, p. 37, Koop)    
    dofs_posterior=T+dofs_prior;
    dofd_posterior=T+dofd_prior;
    
    % Variables for storage
    Sum_supply=zeros(T,1);
    Sum_demand=zeros(T,1);
    Sum_lamda=zeros(T,1);
    
    
    
    %% Burn-off phase
    n_iter=5000;
    
    std_factor=1;
    
    Coeffs_previous=Coeffs_OLS;
    Coeffd_previous=Coeffd_OLS;   
    std_s_previous=std_factor*sqrt(stats_s(4));
    std_d_previous=std_factor*sqrt(stats_d(4));
    
    for in=1:n_iter
        
        lamda=1-(normcdf((RHS_d*Coeffd_previous-RHS_s*Coeffs_previous)./(sqrt(std_s_previous^2+std_d_previous^2)),0,1));
        u=rand(T,1);
        i_d=find(u<lamda);        
        i_s=find(u>=lamda);
        nd=size(i_d,1);
        ns=size(i_s,1);
        yd(i_d)=LHS(i_d);
        ys(i_s)=LHS(i_s);
        mean_s_augm=RHS_s*Coeffs_previous;
        std_s_augm=std_s_previous*ones(T,1);
        mean_d_augm=RHS_d*Coeffd_previous;
        std_d_augm=std_d_previous*ones(T,1);
        
        %Fill augmented data into yd and ys
        for ii=1:nd
            ys(i_d(ii))=randraw_a('normaltrunc',[LHS(i_d(ii)),inf,mean_s_augm(i_d(ii)),std_s_augm(i_d(ii))],1); 
        end;
        for ii=1:ns
            yd(i_s(ii))=randraw_a('normaltrunc',[LHS(i_s(ii)),inf,mean_d_augm(i_s(ii)),std_d_augm(i_s(ii))],1);
        end;
        

        if isempty(find(isnan(yd)))==0
            boomyd
        end;
        if isempty(find(isnan(ys)))==0
            boomys
        end;
        
        % OLS estimation with the augmented data
        % It is used in order to define parameters of the posterior (e.g.
        % mean value of coefficients)
        [Coeffs_OLS,bint_s,r_s,rint_s,stats_s] = regress(ys,RHS_s);
        [Coeffd_OLS,bint_r,r_d,rint_d,stats_d] = regress(yd,RHS_d);
        
        % Calculate OLS variances (Eq. 3.6, p. 36, Koop)        
        ss2_OLS=((ys-RHS_s*Coeffs_OLS)'*(ys-RHS_s*Coeffs_OLS))/dofs_OLS;
        sd2_OLS=((yd-RHS_d*Coeffd_OLS)'*(yd-RHS_d*Coeffd_OLS))/dofd_OLS;
                
        % Calculate posterior variances (Eq. 3.13, p. 37, Koop)
        ss2=(dofs_prior*ss2_prior+dofs_OLS*ss2_OLS+(Coeffs_OLS-Coeffs_prior)'*inv(Var_prior_s+inv(RHS_s'*RHS_s))*(Coeffs_OLS-Coeffs_prior))/dofs_posterior;
        sd2=(dofd_prior*sd2_prior+dofd_OLS*sd2_OLS+(Coeffd_OLS-Coeffd_prior)'*inv(Var_prior_d+inv(RHS_d'*RHS_d))*(Coeffd_OLS-Coeffd_prior))/dofd_posterior;
        
        if isnan(ss2)||isnan(sd2)||isinf(ss2)||isinf(sd2)||ss2>1e03
            boommm
        end;
        
        % Calculate posterior mean values of coefficients (Eq. 3.11, p. 37,
        % Koop)
        Coeffs_posterior=Vs_posterior*(inv(Var_prior_s)*Coeffs_prior+XXs*Coeffs_OLS);
        Coeffd_posterior=Vd_posterior*(inv(Var_prior_d)*Coeffd_prior+XXd*Coeffd_OLS);
        
        % Drawing from the posterior        
            accept = 0;
            
            % Draw from the marginal posterior for h (from a gamma distr) and then find std_s
            h_s=gamrnd(1/sqrt(ss2),dofs_posterior);
            std_s_previous=sqrt(1/h_s);       
            
            % Draw from the conditional posterior (given h) for the coefficients
            % (See (Definition B.26, p. 329, Koop)  and (Eq. 3.9, p. 37, Koop)
            while accept == 0; % rejection bounds on rho
                
                Coeffs_current=Coeffs_posterior+Chol_s_mat*randn(Ns+1,1);
                
                if isnan(Coeffs_current)
                    Coeffs_current                    
                end;
                
                check_supply_restrictions=find((supply_restrictions.*Coeffs_current)<0);
                
                if isempty(check_supply_restrictions)
                    Coeffs_burnoff(in,:)= Coeffs_current';
                    stds_mat_burnoff(in)=std_s_previous;
                    ss2_mat_burnoff(in)=ss2;
                    accept=1;
                end;
            end;
                
            
            accept = 0;
            
            % Draw from the marginal posterior for h (from a gamma distr)
            % and then find std_s
            h_d=gamrnd(1/sqrt(sd2),dofd_posterior);
            std_d_previous=sqrt(1/h_d);       
            
            % Draw from the conditional posterior (given h) for the
            % coefficients
            while accept == 0; % rejection bounds on rho
                
                Coeffd_current=Coeffd_posterior+Chol_d_mat*randn(Nd+1,1);
                
                if isnan(Coeffd_current)
                    Coeffd_current                    
                end;                
                
                check_demand_restrictions=find((demand_restrictions.*Coeffd_current)<0);
                                
                if isempty(check_demand_restrictions)
                    Coeffd_burnoff(in,:)= Coeffd_current';
                    stdd_mat_burnoff(in)=std_d_previous;
                    sd2_mat_burnoff(in)=sd2;
                    accept=1;                
                end;
            end;
                    
            
            % Store data
            ys_mat_burnoff(in,:)=ys';
            yd_mat_burnoff(in,:)=yd'; 
                    
            esupply_mat_burnoff(in,:)=RHS_s*Coeffs_previous;
            edemand_mat_burnoff(in,:)=RHS_d*Coeffd_previous;
            
            lamda_mat_burnoff(in,:)=lamda;           
            
            
            % Update of coefficients
            Coeffs_previous=Coeffs_current;
            Coeffd_previous=Coeffd_current;
            
    end;
    
    
    %% Main loop
    n_iter=20000;

        
    for in=1:n_iter
        
        % Data augmentation
            % lamda is the probability of demand being less than supply,
            % i.e. we are in a demand regime
        lamda=1-(normcdf((RHS_d*Coeffd_previous-RHS_s*Coeffs_previous)./(sqrt(std_s_previous^2+std_d_previous^2)),0,1));
        u=rand(T,1);
        i_d=find(u<lamda);        
        i_s=find(u>=lamda);
        nd=size(i_d,1);
        ns=size(i_s,1);
        yd(i_d)=LHS(i_d);
        ys(i_s)=LHS(i_s);
        mean_s_augm=RHS_s*Coeffs_previous;
        std_s_augm=std_s_previous*ones(T,1);
        mean_d_augm=RHS_d*Coeffd_previous;
        std_d_augm=std_d_previous*ones(T,1);
        
        % Now i have to define yd for the i_s indices and ys for the i_d
        % indices        
        
        %Fill augmented data into yd and ys
        for ii=1:nd
            ys(i_d(ii))=randraw_a('normaltrunc',[LHS(i_d(ii)),inf,mean_s_augm(i_d(ii)),std_s_augm(i_d(ii))],1);            
        end;
        for ii=1:ns
            yd(i_s(ii))=randraw_a('normaltrunc',[LHS(i_s(ii)),inf,mean_d_augm(i_s(ii)),std_d_augm(i_s(ii))],1);
        end;

        
        
        if isempty(find(isnan(yd)))==0
            boomyd
        end;
        if isempty(find(isnan(ys)))==0
            boomys
        end;
        
        
        % 
        % OLS estimation with the augmented data
        % It is used in order to define parameters of the posterior (e.g.
        % mean value of coefficients)
        [Coeffs_OLS,bint_s,r_s,rint_s,stats_s] = regress(ys,RHS_s);
        [Coeffd_OLS,bint_r,r_d,rint_d,stats_d] = regress(yd,RHS_d);

        % Calculate OLS variances (Eq. 3.6, p. 36, Koop)   
        ss2_OLS=((ys-RHS_s*Coeffs_OLS)'*(ys-RHS_s*Coeffs_OLS))/dofs_OLS;
        sd2_OLS=((yd-RHS_d*Coeffd_OLS)'*(yd-RHS_d*Coeffd_OLS))/dofd_OLS;
        
        % Calculate posterior variances (Eq. 3.13, p. 37, Koop)
        ss2=(dofs_prior*ss2_prior+dofs_OLS*ss2_OLS+(Coeffs_OLS-Coeffs_prior)'*inv(Var_prior_s+inv(RHS_s'*RHS_s))*(Coeffs_OLS-Coeffs_prior))/dofs_posterior;
        sd2=(dofd_prior*sd2_prior+dofd_OLS*sd2_OLS+(Coeffd_OLS-Coeffd_prior)'*inv(Var_prior_d+inv(RHS_d'*RHS_d))*(Coeffd_OLS-Coeffd_prior))/dofd_posterior;

        
        % Calculate posterior mean values of coefficients (Eq. 3.11, p. 37,
        % Koop)
        Coeffs_posterior=Vs_posterior*(inv(Var_prior_s)*Coeffs_prior+XXs*Coeffs_OLS);
        Coeffd_posterior=Vd_posterior*(inv(Var_prior_d)*Coeffd_prior+XXd*Coeffd_OLS);
        
        
        % Drawing from the posterior        
            accept = 0;
            
            % Draw from the marginal posterior for h (from a gamma distr)
            % and then find std_s
            h_s=gamrnd(1/sqrt(ss2),dofs_posterior);
            std_s_previous=sqrt(1/h_s);   
            
            % Draw from the conditional posterior (given h) for the
            % coefficients
            while accept == 0; % rejection bounds on rho
                
                Coeffs_current=Coeffs_posterior+Chol_s_mat*randn(Ns+1,1);
                
                if isnan(Coeffs_current)
                    Coeffs_current                    
                end;

                
                check_supply_restrictions=find((supply_restrictions.*Coeffs_current)<0);
                
                if isempty(check_supply_restrictions)
                    Coeffs(in,:)= Coeffs_current';
                    stds_mat(in)=std_s_previous;
                    ss2_mat(in)=ss2;
                    accept=1;
                end;
            end;
                
            
            accept = 0;
            
            % Draw from the marginal posterior for h (from a gamma distr)
            % and then find std_s
            h_d=gamrnd(1/sqrt(sd2),dofd_posterior);
            std_d_previous=sqrt(1/h_d);
            
            % Draw from the conditional posterior (given h) for the
            % coefficients
            while accept == 0; % rejection bounds on rho
                
                Coeffd_current=Coeffd_posterior+Chol_d_mat*randn(Nd+1,1);


                if isnan(Coeffd_current)
                    Coeffd_current                    
                end;                
                
                
                check_demand_restrictions=find((demand_restrictions.*Coeffd_current)<0);
                                
                if isempty(check_demand_restrictions)
                    Coeffd(in,:)= Coeffd_current';
                    stdd_mat(in)=std_d_previous;
                    sd2_mat(in)=sd2;
                    accept=1;                
                end;
            end;
            
        % Storing useful quantities
            
            Sum_supply=Sum_supply+ys;
            Sum_demand=Sum_demand+yd;
            Sum_lamda=Sum_lamda+lamda;
        
            % Store data
            ys_mat(in,:)=ys';
            yd_mat(in,:)=yd'; 
        
            esupply_mat(in,:)=RHS_s*Coeffs_previous;
            edemand_mat(in,:)=RHS_d*Coeffd_previous;
            
            lamda_mat(in,:)=lamda;
            
            i_yd_regime=find(ys>yd);
            ys_rel_yd(in,i_yd_regime)=0;
            i_ys_regime=find(ys<=yd);
            ys_rel_yd(in,i_ys_regime)=1;
            
        % Update of coefficients
            Coeffs_previous=Coeffs_current;
            Coeffd_previous=Coeffd_current;
      
    end;

    
    Mean_demand=Sum_demand/n_iter;
    Mean_supply=Sum_supply/n_iter;
    Mean_lamda=Sum_lamda/n_iter;
    
    for it=1:T
        i_no_outliers_s=find(ys_mat(:,it)<(LHS(it)+4*mean(stds_mat)));
        Mean_supply_2(it)=median(ys_mat(i_no_outliers_s,it));
    end;
    
    for it=1:T
        i_no_outliers_d=find(yd_mat(:,it)<(LHS(it)+4*mean(stdd_mat)));
        Mean_demand_2(it)=median(yd_mat(i_no_outliers_d,it));
    end;
    
    %% Convergence diagnostics
    
    % Diagnostic 1
    % -------------------
    % Diagnostic suggested by Geweke (1992) - check if the estimates based
    % on the first half of the draws is essentially the same as the
    % estimate based on the last half
    
    %vnames_s=strvcat('beta_s1','beta_s2');
    %vnames_s=strvcat('beta_s1','beta_s2','beta_s3');
    %vnames_s=strvcat('beta_s1','beta_s2','beta_s3','beta_s4');
    %vnames_s=strvcat('beta_s1','beta_s2','beta_s3','beta_s4','beta_s5');
    vnames_s=strvcat('beta_s1','beta_s2','beta_s3','beta_s4','beta_s5','beta_s6');
    %vnames_s=strvcat('beta_s1','beta_s2','beta_s3','beta_s4','beta_s5','beta_s6','beta_s7');
    %result_supply=coda(Coeffs,vnames_s);
    %coda(Coeffs,vnames_s);
    
    vnames_d=strvcat('beta_d1','beta_d2','beta_d3');
    %vnames_d=strvcat('beta_d1','beta_d2','beta_d3','beta_d4');
    %vnames_d=strvcat('beta_d1','beta_d2','beta_d3','beta_d4','beta_d5');
    %result_demand=coda(Coeffd,vnames_d);
    %coda(Coeffd,vnames_d);

    %% Results
    
    
        % Set up dates for the plots
    
        years_mat_num=[2003*ones(1,1) 2004*ones(1,12) 2005*ones(1,12) 2006*ones(1,12) 2007*ones(1,12) 2008*ones(1,12) 2009*ones(1,12) 2010*ones(1,12) 2011*ones(1,3)]';
        years_mat=int2str(years_mat_num);
        
        months_mat_num=[12:12 1:12 1:12 1:12 1:12 1:12 1:12 1:12 1:3]';
        months_mat=int2str(months_mat_num);
        
        datess=strcat(months_mat,'-',years_mat);
        
        
        dates_mat=datenum(datess,'mm-yyyy');
        
        
    
        %--------------------------
    
        mult_factor_std=1;
        
        figure;
        for i=1:Ns+1
            mean_temp=mean(Coeffs(:,i));
            std_temp=std(Coeffs(:,i));
            Coeffs(:,i)=(Coeffs(:,i)-mean_temp).*mult_factor_std+mean_temp;            
            subplot(Ns+2,1,(i-1)*1+1); hist(Coeffs(:,i),30);
            ht=title(strcat('Posterior density of b2',num2str(i)));  
            set(ht,'Fontsize',16);   
            haxes=get(gcf,'CurrentAxes');
            set(haxes,'FontSize',16);            
        end;
            subplot(Ns+2,1,(i)*1+1); hist(stds_mat,30);
            ht=title('Posterior density of standard deviation');
            set(ht,'Fontsize',16);   
            haxes=get(gcf,'CurrentAxes');
            set(haxes,'FontSize',16);            
       
        
         figure;
        for i=1:Nd+1
             mean_temp=mean(Coeffd(:,i));
             std_temp=std(Coeffd(:,i));
             Coeffd(:,i)=(Coeffd(:,i)-mean_temp).*mult_factor_std+mean_temp;
            subplot(Nd+2,1,(i-1)*1+1); hist(Coeffd(:,i),30);
            ht=title(strcat('Posterior density of b1',num2str(i)));  
            set(ht,'Fontsize',16);   
            haxes=get(gcf,'CurrentAxes');
            set(haxes,'FontSize',16);            
        end;
            subplot(Nd+2,1,(i)*1+1); hist(stdd_mat,30);
            ht=title('Posterior density of standard deviation');
            set(ht,'Fontsize',16);   
            haxes=get(gcf,'CurrentAxes');
            set(haxes,'FontSize',16);                    
        
      
% Means and quartiles

    ysmean=mean(ys_mat,1); 
    ydmean=mean(yd_mat,1);
    figure;
    plot(dates_mat,ysmean,'r','LineWidth',2);
    hold on;
    ysmean75 = prctile(ys_mat,90,1);
    hs75=plot(dates_mat,ysmean75,'--r','LineWidth',1.5);
    hAnnotations75 = get(hs75,'Annotation');
    hLegendEntrys75 = get(hAnnotations75','LegendInformation');
    set(hLegendEntrys75,'IconDisplayStyle','off') 
    hold on;
    ysmean25 = prctile(ys_mat,10,1);
    hs25=plot(dates_mat,ysmean25','-.r','LineWidth',1.5);
    hAnnotations25 = get(hs25,'Annotation');
    hLegendEntrys25 = get(hAnnotations25','LegendInformation');
    set(hLegendEntrys25,'IconDisplayStyle','off') 
    hold on;
    plot(dates_mat,ydmean,'b','LineWidth',2);
    hold on;
    ydmean75 = prctile(yd_mat,90,1);
    hd75=plot(dates_mat,ydmean75,'--b','LineWidth',1.5);
    hAnnotationd75 = get(hd75,'Annotation');
    hLegendEntryd75 = get(hAnnotationd75','LegendInformation');
    set(hLegendEntryd75,'IconDisplayStyle','off') 
    hold on;
    ydmean25 = prctile(yd_mat,10,1);
    hd25=plot(dates_mat,ydmean25','-.b','LineWidth',1.5);
    hAnnotationd25 = get(hd25,'Annotation');
    hLegendEntryd25 = get(hAnnotationd25','LegendInformation');
    set(hLegendEntryd25,'IconDisplayStyle','off')       
    hold on;
    
    plot(dates_mat,LHS,'k','LineWidth',1);     
     
    set(gca,'XTick',731916:180:734563)
    xlim([731916 734563]);
    datetick('x','mm/yy','keepticks','keeplimits');        
    legend('Supply','Demand','Observed credit');
    haxes=get(gcf,'CurrentAxes');
    set(haxes,'FontSize',18);  
    
    ymin=11;
    ymax=13;
    ylim([ymin ymax]);

    n1 = datenum('04/10', 'mm/yy');
    n2 = datenum('04/10', 'mm/yy');
    l1=line([n1 n1], [ymin ymax]);
    set(l1,'Color','k');
    set(l1,'LineWidth',1.5);
    t1=text(n1,11.3,'FSM');
    set(t1,'FontSize',16);
    
    n1 = datenum('09/08', 'mm/yy');
    n2 = datenum('09/08', 'mm/yy');
    l1=line([n1 n1], [ymin ymax]);
    set(l1,'Color','k');
    set(l1,'LineWidth',1.5);
    t1=text(n1,11.3,'Lehman collapse'); 
    set(t1,'FontSize',16);
    
    %------------------------------------
    
    
    figure;
    plot(dates_mat,ysmean-LHS','r','LineWidth',2);
    hold on;    
    yss75=plot(dates_mat,ysmean75-LHS','--r','LineWidth',1);
    hhAnnotations75 = get(yss75,'Annotation');
    hhLegendEntrys75 = get(hhAnnotations75','LegendInformation');
    set(hhLegendEntrys75,'IconDisplayStyle','off') 
    hold on;
    yss25=plot(dates_mat,ysmean25-LHS','--r','LineWidth',1);
    hhAnnotations25 = get(yss25,'Annotation');
    hhLegendEntrys25 = get(hhAnnotations25','LegendInformation');
    set(hhLegendEntrys25,'IconDisplayStyle','off')     
    hold on;
    plot(dates_mat,ydmean-LHS','b','LineWidth',2);
    hold on;    
    ydd25=plot(dates_mat,ydmean25-LHS','--b','LineWidth',1);
    hhhAnnotations25 = get(ydd25,'Annotation');
    hhhLegendEntrys25 = get(hhhAnnotations25','LegendInformation');
    set(hhhLegendEntrys25,'IconDisplayStyle','off')         
    hold on;
    ydd75=plot(dates_mat,ydmean75-LHS','--b','LineWidth',1);
    hhhAnnotations75 = get(ydd75,'Annotation');
    hhhLegendEntrys75 = get(hhhAnnotations75','LegendInformation');
    set(hhhLegendEntrys75,'IconDisplayStyle','off')         
    hold on;
    h0=plot(dates_mat,LHS-LHS,'k','LineWidth',1);  
    hAnnotation0 = get(h0,'Annotation');
    hLegendEntry0 = get(hAnnotation0','LegendInformation');
    set(hLegendEntry0,'IconDisplayStyle','off')       
    set(gca,'XTick',731916:180:734563)
    xlim([731916 734563]);
    datetick('x','mm/yy','keepticks','keeplimits');        
    legend('Excess supply','Excess demand');
    haxes=get(gcf,'CurrentAxes');
    set(haxes,'FontSize',18);
        
    ymin=-0.2;
    ymax=0.5;
    ylim([ymin ymax]);
    
    n1 = datenum('04/10', 'mm/yy');
    n2 = datenum('04/10', 'mm/yy');
    l1=line([n1 n1], [ymin ymax]);
    set(l1,'Color','k');
    set(l1,'LineWidth',1.5);
    t1=text(n1,0.05,'FSM');
    set(t1,'FontSize',16);
    
    n1 = datenum('09/08', 'mm/yy');
    n2 = datenum('09/08', 'mm/yy');
    l1=line([n1 n1], [ymin ymax]);
    set(l1,'Color','k');
    set(l1,'LineWidth',1.5);
    t1=text(n1,0.05,'Lehman collapse'); 
    set(t1,'FontSize',16);
    
                
    lamda_mean=mean(lamda_mat,1);
    lamda_max=prctile(lamda_mat,75,1); 
    lamda_min=prctile(lamda_mat,25,1); 
    lamda_final=lamda_mean; 
    figure;
    plot(dates_mat,lamda_mean,'r','LineWidth',2);   
    hold on;    
    plot(dates_mat,lamda_max,'--r','LineWidth',1.5);
    hold on;    
    plot(dates_mat,lamda_min,'--r','LineWidth',1.5);
    set(gca,'XTick',731916:180:734563)
    xlim([731916 734563]);
    datetick('x','mm/yy','keepticks','keeplimits');        
    legend('Mean probability of being in demand regime (bounded by the first and third quartile)');
    haxes=get(gcf,'CurrentAxes');
    set(haxes,'FontSize',18);
    
    n1 = datenum('04/10', 'mm/yy');
    n2 = datenum('04/10', 'mm/yy');
    l1=line([n1 n1], [0 1]);
    set(l1,'Color','k');
    set(l1,'LineWidth',1.5);
    t1=text(n1,0.5,'FSM'); 
    set(t1,'FontSize',16);
    
    
    n1 = datenum('09/08', 'mm/yy');
    n2 = datenum('09/08', 'mm/yy');
    l1=line([n1 n1], [0 1]);
    set(l1,'Color','k');
    set(l1,'LineWidth',1.5);
    t1=text(n1,0.5,'Lehman collapse'); 
    set(t1,'FontSize',16);
    
    
    
    
    
    
    %-----------------------------
    
     prob_ysyd=1-mean(ys_rel_yd,1);
    
    prob_ysyd_final=prob_ysyd; 
    figure;
    plot(dates_mat,prob_ysyd_final,'r');   
    set(gca,'XTick',731916:180:734563)
    xlim([731916 734563]);
    datetick('x','mm/yy','keepticks','keeplimits');        
    legend('Mean probability of being in demand regime');
    haxes=get(gcf,'CurrentAxes');
    set(haxes,'FontSize',16);
    
    n1 = datenum('04/10', 'mm/yy');
    n2 = datenum('04/10', 'mm/yy');
    l1=line([n1 n1], [0 1]);
    set(l1,'Color','k');
    set(l1,'LineWidth',1.5);
    t1=text(n1,0.8*ymax/(ymax-ymin),'FSM'); 
    set(t1,'FontSize',14);
    
    
    n1 = datenum('09/08', 'mm/yy');
    n2 = datenum('09/08', 'mm/yy');
    l1=line([n1 n1], [0 1]);
    set(l1,'Color','k');
    set(l1,'LineWidth',1.5);
    t1=text(n1,0.8*ymax/(ymax-ymin),'Lehman collapse'); 
    set(t1,'FontSize',14);
    
    
    %% OTHER PLOTS
    %----------
    
    
    
    %   -----  plot ecb dependence
    figure;
    
    plot(dates_mat,ecb_dependence,'r');   
    set(gca,'XTick',731916:180:734563)
    xlim([731916 734563]);
    datetick('x','mm/yy','keepticks','keeplimits');     
    title('ECB dependence');
    haxes=get(gcf,'CurrentAxes');
    set(haxes,'FontSize',16);
    
    n1 = datenum('04/10', 'mm/yy');
    n2 = datenum('04/10', 'mm/yy');
    l1=line([n1 n1], [0 1]);
    set(l1,'Color','k');
    set(l1,'LineWidth',1.5);
    t1=text(n1,0.2,'FSM'); 
    set(t1,'FontSize',12);
    
    
    n1 = datenum('09/08', 'mm/yy');
    n2 = datenum('09/08', 'mm/yy');
    l1=line([n1 n1], [0 1]);
    set(l1,'Color','k');
    set(l1,'LineWidth',1.5);
    t1=text(n1,0.2,'Lehman collapse'); 
    set(t1,'FontSize',12);
    
    
