% dols.m: Uses dynamic OLS to estimate the cointegrating relationships
% considered in my paper, "On the Welfare Cost of Inflation and the Recent
% Behavior of Money Demand,"  so as to produce the numbers in Tables 4 and 
% 5.
%
%   Takes as input the data from quarterly.dat, containing post-1980
%   quarterly series. Column 1 is real money balances, column 2 is real
%   GDP, column 3 is the money-income ratio, and column 4 is the nominal
%   interest rate.
%
%   Produces as output a set of matrices betamatx, sematx, and tmatx or
%   chi2matx containing the regression cofficients, standard errors, and t
%   or chi^2 statistics for testing the null hypotheses of a negative
%   interest elasticity or a unitary income elasticity of money demand,
%   where x = lmr for the regression linking the log of the money-income
%   ratio to the level of the nominal interest rate and x = lMPlYpr for the
%   regression linking the log of real balances to the log of real GDP and
%   the level of the nominal interest rate. In each case, these results for
%   presented for all combinations of values of p ranging between 1 and 4
%   8 and q ranging between 0 and 8, where p is the number of leads and lags
%   included in the DOLS regression and q is the Newey-West lag truncation
%   parameter.
%
%   Also produces as output matrices outmatx for the same x's indicated
%   above, which reproduce the numbers shown in Tables 4 and 5 of the 
%   paper. For the regressions involving the money-income ratio, 
%   columns 1 and 2 keep track of p and q, columns 3 and 4 display the
%   regression coefficients on the constant and the interest rate term,
%   columns 5 and 6 display the corresponding standard errors, and column 7
%   displays the t-statistic for testing for a negative interest elasticity
%   of money demand. For the regressions involving real money balances, 
%   columns 1 and 2 keep track of p and q, columns 3, 4, and 5 display the
%   regression coefficients on the constant, real GDP, and the interest
%   rate term, columns 6, 7, and 8 display the corresponding standard
%   errors, and column 9 displays the chi^2 statistic for testing for a
%   unitary income elasticity of money demand.
%
%   Written for MATLAB Version 7.6 (2008a) by
%
%     Peter N. Ireland
%     Boston College
%     Department of Economics
%     140 Commonwealth Avenue
%     Chestnut Hill, MA 02467 USA
%     irelandp@bc.edu
%
%   COPYRIGHT (c) 2008 BY PETER N. IRELAND.  REDISTRIBUTION IS PERMITTED FOR
%   EDUCATIONAL AND RESEARCH PURPOSES, SO LONG AS NO CHANGES ARE MADE. ALL
%   COPIES MUST BE PROVIDED FREE OF CHARGE AND MUST INCLUDE THIS COPYRIGHT
%   NOTICE.

% estimate equations with lm and r

  load quarterly.dat;
  
  yt = log(quarterly(:,3));
  
  T = length(yt);
  
  xt = [ quarterly(:,4) ];
  
  dxt = xt(2:T,:) - xt(1:T-1,:);
  
  betamatlmr = zeros(36,11);
  
  sematlmr = zeros(36,11);
  
  tmatlmr = zeros(36,1);
  
  for p = 1:4
      
      bigyt = yt(2+p:T-p);
      
      bigxt = zeros(T-1-2*p,3+2*p);
      
      bigxt(:,1:3) = [ ones(T-1-2*p,1) xt(2+p:T-p,:) dxt(1+p:T-1-p,:) ];
      
      if p >= 1
          
          for pp = 1:p
              
              bigxt(:,2*pp+2:2*pp+3) = [ dxt(1+p-pp:T-1-p-pp,:) dxt(1+p+pp:T-1-p+pp,:) ];
              
          end
          
      end
      
      bigbeta = inv(bigxt'*bigxt)*bigxt'*bigyt;
      
      ut = bigyt - bigxt*bigbeta;
      
      s2 = (1/(T-4-4*p))*ut'*ut;
      
      gamo = (1/(T-1-2*p))*ut'*ut;

      gamvec = zeros(8,1);

      for j = 1:8

        gamvec(j) = (1/T)*ut(j+1:T-1-2*p)'*ut(1:T-1-2*p-j);

      end
      
      betamatlmr(9*p-8,1:2*p+3) = bigbeta';
  
      sematlmr(9*p-8,1:2*p+3) = sqrt(s2*diag(inv(bigxt'*bigxt)));
      
      tmatlmr(9*p-8) = betamatlmr(9*p-8,2)/sematlmr(9*p-8,2);
      
      for q = 1:8
          
         lam2  = gamo;
         
         for j = 1:q
             
             lam2 = lam2 + 2*(1-j/(q+1))*gamvec(j);
             
         end
         
         betamatlmr(9*p-8+q,1:2*p+3) = bigbeta';
  
         sematlmr(9*p-8+q,1:2*p+3) = sqrt(lam2*diag(inv(bigxt'*bigxt)));
         
         tmatlmr(9*p-8+q) = betamatlmr(9*p-8+q,2)/sematlmr(9*p-8+q,2);
          
      end
      
  end

  outmatlmr = zeros(36,7);
  
  for p = 1:4
        
       for q = 0:8
            
            outmatlmr(9*p-8+q,:) = [ p q betamatlmr(9*p-8+q,1:2) sematlmr(9*p-8+q,1:2) tmatlmr(9*p-8+q) ];
            
        end
        
  end
  
  outmatlmr

% estimate equations with lMP, lYP, and r
  
  yt = log(quarterly(:,1));
  
  T = length(yt);
  
  xt = [ log(quarterly(:,2)) quarterly(:,4) ];
  
  dxt = xt(2:T,:) - xt(1:T-1,:);
  
  betamatlMPlYPr = zeros(36,21);
  
  sematlMPlYPr = zeros(36,21);
  
  chi2matlMPlYPr = zeros(36,1);
  
  for p = 1:4
      
      bigyt = yt(2+p:T-p);
      
      bigxt = zeros(T-1-2*p,5+4*p);
      
      bigxt(:,1:5) = [ ones(T-1-2*p,1) xt(2+p:T-p,:) dxt(1+p:T-1-p,:) ];
      
      if p >= 1
          
          for pp = 1:p
              
              bigxt(:,4*pp+2:4*pp+5) = [ dxt(1+p-pp:T-1-p-pp,:) dxt(1+p+pp:T-1-p+pp,:) ];
              
          end
          
      end
      
      bigbeta = inv(bigxt'*bigxt)*bigxt'*bigyt;
      
      ut = bigyt - bigxt*bigbeta;
      
      s2 = (1/(T-6-6*p))*ut'*ut;
      
      gamo = (1/(T-1-2*p))*ut'*ut;

      gamvec = zeros(8,1);

      for j = 1:8

        gamvec(j) = (1/T)*ut(j+1:T-1-2*p)'*ut(1:T-1-2*p-j);

      end
      
      betamatlMPlYPr(9*p-8,1:4*p+5) = bigbeta';
  
      sematlMPlYPr(9*p-8,1:4*p+5) = sqrt(s2*diag(inv(bigxt'*bigxt)));
  
      chi2matlMPlYPr(9*p-8) = ((betamatlMPlYPr(9*p-8,2)-1)/sematlMPlYPr(9*p-8,2))^2;
      
      for q = 1:8
          
         lam2  = gamo;
         
         for j = 1:q
             
             lam2 = lam2 + 2*(1-j/(q+1))*gamvec(j);
             
         end
         
         betamatlMPlYPr(9*p-8+q,1:4*p+5) = bigbeta';
  
         sematlMPlYPr(9*p-8+q,1:4*p+5) = sqrt(lam2*diag(inv(bigxt'*bigxt)));
  
         chi2matlMPlYPr(9*p-8+q) = ((betamatlMPlYPr(9*p-8+q,2)-1)/sematlMPlYPr(9*p-8+q,2))^2;
          
      end
      
  end

  outmatlMPlYPr = zeros(36,9);
  
  for p = 1:4
        
       for q = 0:8
            
            outmatlMPlYPr(9*p-8+q,:) = [ p q betamatlMPlYPr(9*p-8+q,1:3) sematlMPlYPr(9*p-8+q,1:3) chi2matlMPlYPr(9*p-8+q) ];
            
        end
        
  end
  
  outmatlMPlYPr