% philou.m: Performs the Phillips-Ouliaris tests for cointegration to the
% 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 2 and 3.
%
%   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 postatx, containing
%   values of the Phillips-Ouliaris statistic for values of the Newey-West
%   lag trunction parameter q between 0 and 8, where x = lmlr for the
%   regression linking the log of the money-income ratio to the log of the
%   nominal interest rate, x = lmr for the regression linking the log of
%   the money-income ratio to the level of the nominal interest rate, and 
%   x = LMPLYPlr for the regression linking the log of real balances to the
%   log of real GDP and the log of the nominal interest rate. A constant is
%   included in each regression.
%
%   Also produces as output matrices outmatx for the same x's indicated
%   above, which reproduce the numbers shown in Tables 2 and 3 of the
%   paper. For x = lmlr and x = lmr, the top row of columns 1 through 3
%   displays the regression coefficients and the serial correlation
%   parameter for the regression residuals, and the rows of columns 4 and 5
%   display the Phillips-Ouliaris statistic (column 5) for each value of q
%   (column 4). For x = lMPlYPlr, the top row of columns 1 through 4
%   displays the regression coefficients and the serial correlation
%   parameter for the regression residuals, and the rows of columns 5 and 6
%   display the Phillips-Ouliaris statistic (column 6) for each value of q
%   (column 5).
%
%   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 cointegration between lm and lr

  load quarterly.dat;
  
  yt = log(quarterly(:,3));
  
  T = length(yt);
  
  xt = [ ones(T,1) log(quarterly(:,4)) ];

% perform conintegrating regression and obtain residuals

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

  ut = yt - xt*betvec;

% perform residual regression and compute preliminary statistics

  yt = ut(2:T);

  xt = ut(1:T-1);
  
  rho = inv(xt'*xt)*xt'*yt;

  et = yt - xt*rho;

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

  s = sqrt(s2);

  sigr = s*sqrt(inv(xt'*xt));

  tstat = (rho-1)/sigr;

% compute Newey-West variance

  gamo = (1/(T-1))*et'*et;

  gamvec = zeros(50,1);

  for j = 1:51

    gamet = et(j+1:T-1);

    gamej = et(1:T-j-1);

    gamvec(j) = (1/(T-1))*gamet'*gamej;

  end

% compute Phillips-Ouliaris statistic for various values of q

  postatlmlr = zeros(9,1);

  postatlmlr(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);

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

  end

% report results

  qvec  = 0:8;
  
  qvec = qvec';

  outmatlmlr = [ betvec' rho qvec(1) postatlmlr(1) ; ...
                 zeros(8,3) qvec(2:9) postatlmlr(2:9) ];
             
  outmatlmlr
  
% test for cointegration between lm and r
  
  yt = log(quarterly(:,3));
  
  T = length(yt);
  
  xt = [ ones(T,1) quarterly(:,4) ];

% perform conintegrating regression and obtain residuals

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

  ut = yt - xt*betvec;

% perform residual regression and compute preliminary statistics

  yt = ut(2:T);

  xt = ut(1:T-1);
  
  rho = inv(xt'*xt)*xt'*yt;

  et = yt - xt*rho;

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

  s = sqrt(s2);

  sigr = s*sqrt(inv(xt'*xt));

  tstat = (rho-1)/sigr;

% compute Newey-West variance

  gamo = (1/(T-1))*et'*et;

  gamvec = zeros(50,1);

  for j = 1:51

    gamet = et(j+1:T-1);

    gamej = et(1:T-j-1);

    gamvec(j) = (1/(T-1))*gamet'*gamej;

  end

% compute Phillips-Ouliaris statistic for various values of q

  postatlmr = zeros(9,1);

  postatlmr(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);

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

  end

% report results

  outmatlmr = [ betvec' rho qvec(1) postatlmr(1) ; ...
                zeros(8,3) qvec(2:9) postatlmr(2:9) ];
             
  outmatlmr
  
% test for cointegration between lMP lYP and lr

  load quarterly.dat;
  
  yt = log(quarterly(:,1));
  
  T = length(yt);
  
  xt = [ ones(T,1) log(quarterly(:,2)) log(quarterly(:,4)) ];

% perform conintegrating regression and obtain residuals

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

  ut = yt - xt*betvec;

% perform residual regression and compute preliminary statistics

  yt = ut(2:T);

  xt = ut(1:T-1);
  
  rho = inv(xt'*xt)*xt'*yt;

  et = yt - xt*rho;

  s2 = et'*et/(T-3);

  s = sqrt(s2);

  sigr = s*sqrt(inv(xt'*xt));

  tstat = (rho-1)/sigr;

% compute Newey-West variance

  gamo = (1/(T-1))*et'*et;

  gamvec = zeros(50,1);

  for j = 1:51

    gamet = et(j+1:T-1);

    gamej = et(1:T-j-1);

    gamvec(j) = (1/(T-1))*gamet'*gamej;

  end

% compute Phillips-Ouliaris statistic for various values of q

  postatlMPlYPlr = zeros(9,1);

  postatlMPlYPlr(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);

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

  end

% report results

  qvec  = 0:8;
  
  qvec = qvec';

  outmatlMPlYPlr = [ betvec' rho qvec(1) postatlMPlYPlr(1) ; ...
                    zeros(8,4) qvec(2:9) postatlMPlYPlr(2:9) ];
             
  outmatlMPlYPlr
  