/* MSEDIF3      7/21/95 mkh */

/* A program to compute real exchange rate variance decompositions. */
/* Run MSEPLOT right after running this program to graph the results */

 t1=1;    t2=408;           @ whole sample:  1/62-3/95, T=399 @
@  t1=1;    t2=93;  @           @ fixed rate sample:  1/62-9/69, T=93 @
@ t1=133;  t2=399;  @         @ flexible rate sample:  1/73-3/95, T=267 @

load lcpi;
lcpi = lcpi[t1:t2,.];
pcan = lcpi[.,1:2];
pfra = lcpi[.,3:4];
pger = lcpi[.,5:6];
pita = lcpi[.,7:8];
pjap = lcpi[.,9:10];
pusa = lcpi[.,11:12];

/* The price data is now in matrices whose rows are the observations
   at each date.  The first column is the log of the price of traded goods,
   the second column is the log of the price of non-traded goods. */

load weights;
weight = weights;

/* weight is a 6x2 matrix that contains the weights of the components
   of the CPI.  Each row contains the weights for a different country.
   Row 1 is Canada; row 2 is France; row 3 is Germany; row 4 is Italy;
   row 5 is Japan; and, row 6 is the U.S.  The columns contain the weights
   for the two sub-categories of data from the OECD database on
   Datastream.  Column 1 is traded goods and column 2 is non-traded goods.
   The weights were obtained by regressing the log-difference of the CPI
   on the log-differences of the four components over the entire sample */

load lxrate3;
exrat = lxrate3[t1:t2,.];

/*  exrat is a nx5 matrix that contains the time series on exchange rates.
    Each column represents a different dollar price of foreign currency.
    Column 1 is the dollar price of Canadian dollars; column 2 is dollars
    per French franc; column 3 is dollars per German mark; column 4 is
    dollars per Italian lire; column 5 is dollars per Japanese yen.  The
    data are in logs. */

/* Compute x, y and q for each of the 5 countries relative to the U.S. */

yc = weight[6,2]*(pusa[.,2]-pusa[.,1]) - weight[1,2]*(pcan[.,2]-pcan[.,1]);
xc = pusa[.,1] - exrat[.,1] - pcan[.,1];
yf = weight[6,2]*(pusa[.,2]-pusa[.,1]) - weight[2,2]*(pfra[.,2]-pfra[.,1]);
xf = pusa[.,1] - exrat[.,2] - pfra[.,1];
yg = weight[6,2]*(pusa[.,2]-pusa[.,1]) - weight[3,2]*(pger[.,2]-pger[.,1]);
xg = pusa[.,1] - exrat[.,3] - pger[.,1];
yi = weight[6,2]*(pusa[.,2]-pusa[.,1]) - weight[4,2]*(pita[.,2]-pita[.,1]);
xi = pusa[.,1] - exrat[.,4] - pita[.,1];
yj = weight[6,2]*(pusa[.,2]-pusa[.,1]) - weight[5,2]*(pjap[.,2]-pjap[.,1]);
xj = pusa[.,1] - exrat[.,5] - pjap[.,1];

T  = t2 - t1 + 1;            @ number of observations in sample @
tt=seqa(1,1,T);

mc = zeros(T-12,12);         @ Initialize results matrix  for Canada  @
mf = zeros(T-12,12);         @ Initialize results matrix  for France  @
mg = zeros(T-12,12);         @ Initialize results matrix  for Germany @
mi = zeros(T-12,12);         @ Initialize results matrix  for Italy   @
mj = zeros(T-12,12);         @ Initialize results matrix  for Japan   @

/************** Start Procedure *******************************************/
/*
   m = ms(x,y)          Procedure to compute MSE components of x and y

   Inputs       x       Real rate of exchange of US vs foreign traded goods
                y       Weighted difference of US vs foreign relative prices
                        for non-traded goods in terms of traded goods

   Output       m       (T-12)x12 matrix of MSE components for a given
                         foreign country relative to the US
*/

proc (1) = ms(x,y);

local m,n,k,q,qbar,xbar,ybar,qnbar,xnbar,ynbar,var,trend,mse;

q = x+y;

qbar = (1/(T-1))*sumc(q-(q[1]|q[1:T-1]));     @ E(q[t+1]-q[t]) @
xbar = (1/(T-1))*sumc(x-(x[1]|x[1:T-1]));     @ E(x[t+1]-x[t]) @
ybar = (1/(T-1))*sumc(y-(y[1]|y[1:T-1]));     @ E(y[t+1]-y[t]) @

m = zeros(T-12,12);

n = 1;
do until n > T-12;

 k = T/((T-n-1)*(T-n));       @ Factor to correct for small sample bias @

 qnbar = zeros(n,1)|(n*qbar*ones(T-n,1));  @ Tx1 vector with n zeros at top @
 xnbar = zeros(n,1)|(n*xbar*ones(T-n,1));  @ Tx1 vector with n zeros at top @
 ynbar = zeros(n,1)|(n*ybar*ones(T-n,1));  @ Tx1 vector with n zeros at top @

 m[n,2]  = k*moment(q-(q[1:n]|q[1:T-n])-qnbar,0);
 m[n,3]  = (n*qbar)^2;
 m[n,1]  = m[n,2]+m[n,3];
 m[n,5]  = k*moment(x-(x[1:n]|x[1:T-n])-xnbar,0);
 m[n,6]  = (n*xbar)^2;
 m[n,4]  = m[n,5]+m[n,6];
 m[n,8]  = k*moment(y-(y[1:n]|y[1:T-n])-ynbar,0);
 m[n,9]  = (n*ybar)^2;
 m[n,7]  = m[n,8]+m[n,9];
 m[n,11] = k*(x-(x[1:n]|x[1:T-n])-xnbar)'(y-(y[1:n]|y[1:T-n])-ynbar);
 m[n,12] = (n*xbar)*(n*ybar);
 m[n,10] = m[n,11]+m[n,12];

 n = n+1;

endo;

retp(m);
endp;

/************** End Procedure *********************************************/

mc = ms(xc,yc);   @ Call the 'ms' procedure for Can, Fra, Ger, Ita and Jap @
mf = ms(xf,yf);
mg = ms(xg,yg);
mi = ms(xi,yi);
mj = ms(xj,yj);

/* Now run MSEPLOT to generate variance ratio graphs */
/* Make sure the sample range settings are consistent between programs */