% philper.m: Performs the Phillips-Perron unit root tests for the series
% considered in my paper, "On the Welfare Cost of Inflation and the Recent
% Behavior of Money Demand," so as to produce the numbers in Table 1.
%
%   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 9x1 vectors named ppstatx, containg values
%   of the Phillips-Perron statistic for values of the Newey-West lag
%   truncation parameter q between 0 and 8, where x = lm for the log of
%   the money-income ratio, x = lr for the log of the nominal interest
%   rate, and x = r for the level of the nominal interest rate.
%
%   Also produces as output matrices outmatx for the same x's indicated
%   above, which reproduce the numbers shown in Table 1 of the paper: the
%   top row of columns 1 and 2 display the regression coefficients and
%   the rows of columns 3 and 4 display the Phillips-Perron statistic
%   (column 4) for each value of q (column 3).
%
%   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.
  
% test for log money-income ratio

  load quarterly.dat;
 
  yt = log(quarterly(:,3));
  
  T = length(yt) - 1;
  
% set up matrices for regression

  xt = [ ones(T,1) yt(1:T) ];

  yt = yt(2:T+1);

% perform regression and compute preliminary statistics

  rhovec = inv(xt'*xt)*xt'*yt;

  ut = yt - xt*rhovec;

  s2 = ut'*ut/(T-2);

  s = sqrt(s2);

  sigrvec = s2*inv(xt'*xt);

  rho = rhovec(2);

  sigr = sqrt(sigrvec(2,2));

  tstat = (rho-1)/sigr;

% compute Newey-West variance

  gamo = (1/T)*ut'*ut;

  gamvec = zeros(50,1);

  for j = 1:51

    gamu = ut(j+1:T);

    gamuj = ut(1:T-j);

    gamvec(j) = (1/T)*gamu'*gamuj;

  end

% compute Phillips-Perron statistic for various values of q

  ppstatlm = zeros(9,1);

  ppstatlm(1) = tstat;

  for q = 1:8

    lam2 = gamo;

    for j = 1:q

      lam2 = lam2 + 2*(1-j/(q+1))*gamvec(j);

    end

    lam = sqrt(lam2);

    ppstatlm(q+1) = (sqrt(gamo/lam2))*tstat ...
                      - (1/2)*((lam2-gamo)/lam)*(T*sigr/s);

  end
  
% report results

  qvec = 0:8;
   
  qvec = qvec';
  
  outmatlm = [ rhovec' qvec(1) ppstatlm(1) ; ...
               zeros(8,2) qvec(2:9) ppstatlm(2:9) ];
            
  outmatlm
  
% test for log nominal interest rate
 
  yt = log(quarterly(:,4));
  
% set up matrices for regression

  xt = [ ones(T,1) yt(1:T) ];

  yt = yt(2:T+1);

% perform regression and compute preliminary statistics

  rhovec = inv(xt'*xt)*xt'*yt;

  ut = yt - xt*rhovec;

  s2 = ut'*ut/(T-2);

  s = sqrt(s2);

  sigrvec = s2*inv(xt'*xt);

  rho = rhovec(2);

  sigr = sqrt(sigrvec(2,2));

  tstat = (rho-1)/sigr;

% compute Newey-West variance

  gamo = (1/T)*ut'*ut;

  gamvec = zeros(50,1);

  for j = 1:51

    gamu = ut(j+1:T);

    gamuj = ut(1:T-j);

    gamvec(j) = (1/T)*gamu'*gamuj;

  end

% compute Phillips-Perron statistic for various values of q

  ppstatlr = zeros(9,1);

  ppstatlr(1) = tstat;

  for q = 1:8

    lam2 = gamo;

    for j = 1:q

      lam2 = lam2 + 2*(1-j/(q+1))*gamvec(j);

    end

    lam = sqrt(lam2);

    ppstatlr(q+1) = (sqrt(gamo/lam2))*tstat ...
                      - (1/2)*((lam2-gamo)/lam)*(T*sigr/s);

  end
  
% report results
  
  outmatlr = [ rhovec' qvec(1) ppstatlr(1) ; ...
               zeros(8,2) qvec(2:9) ppstatlr(2:9) ];
            
  outmatlr
  
% test for level of nominal interest rate
 
  yt = quarterly(:,4);
  
% set up matrices for regression

  xt = [ ones(T,1) yt(1:T) ];

  yt = yt(2:T+1);

% perform regression and compute preliminary statistics

  rhovec = inv(xt'*xt)*xt'*yt;

  ut = yt - xt*rhovec;

  s2 = ut'*ut/(T-2);

  s = sqrt(s2);

  sigrvec = s2*inv(xt'*xt);

  rho = rhovec(2);

  sigr = sqrt(sigrvec(2,2));

  tstat = (rho-1)/sigr;

% compute Newey-West variance

  gamo = (1/T)*ut'*ut;

  gamvec = zeros(50,1);

  for j = 1:51

    gamu = ut(j+1:T);

    gamuj = ut(1:T-j);

    gamvec(j) = (1/T)*gamu'*gamuj;

  end

% compute Phillips-Perron statistic for various values of q

  ppstatr = zeros(9,1);

  ppstatr(1) = tstat;

  for q = 1:8

    lam2 = gamo;

    for j = 1:q

      lam2 = lam2 + 2*(1-j/(q+1))*gamvec(j);

    end

    lam = sqrt(lam2);

    ppstatr(q+1) = (sqrt(gamo/lam2))*tstat ...
                     - (1/2)*((lam2-gamo)/lam)*(T*sigr/s);

  end
  
% report results
  
  outmatr = [ rhovec' qvec(1) ppstatr(1) ; ...
              zeros(8,2) qvec(2:9) ppstatr(2:9) ];
            
  outmatr