options crt;  
name hptrend 'Hodrick-Prescott Trend Decomposition' ;
?
?  minimize  sum(t=1,T) [y(t) - tau(t)]**2  +
?            sum(t=2,T-1)  lambda*[{tau(t+1)-tau(t)} - {tau(t)-tau(t-1)}]**2
?
?  First order condition for tau(1) to tau(3) [t1 to t3] is:
?   y1 = t1 + lambda*(   t1 - 2*t2 +   t3)
?   y2 = t2 + lambda*(-2*t1 + 5*t2 - 4*t3 +   t4)
?   y3 = t3 + lambda*(   t1 - 4*t2 + 6*t3 - 4*t4 + t5)
?   Y = W*TAU, so TAU = W"Y.  W are the coefficients on t1,t2, etc.
?
?  TSP 4.2B (February 1994)
?  Author: Clint Cummins 
?  using Martin Brookes' matrix method  (slightly revised).
?
load;               ? Get the data.
lrgnp = log(rgnp);
set lambda = 1600;  ? from Finn E. Kydland, Edward C. Prescott
                    ? Business Cycles p.9 (complete ref?)
                    ? N.B. for FREQ M, use lambda >= 14400
                    ?      for FREQ A, use lambda = 100
?
?  Call the HPTREND proc (defined below).
?
HPtrend lrgnp lambda tau;
?
?   plot(dev=lj3,file='hptrend.plt',
?   title='Hodrick-Prescott trend, log(real GNP)') lrgnp tau;
? 
print tau;
?
? =======================================================================
?
proc HPtrend y lambda tau;
  set n = @nob;
  mform(type=sym,nrow=n) w;
?
? Fill in first 2 rows of W by hand:
?
  set w(1) = 1+lambda;
  set w(2) = -2*lambda;
  set ij = n+2;
  set w(ij) = 1+5*lambda;
?
? Fill in most rows of W with a loop:
?
  do i=3,n;
    do imj=0,2;                           ? imj = i-j;  j = i-imj;
      set ij = i + (i-imj-1)*n;
      set w(ij) = (imj=2)*lambda - (imj=1)*4*lambda 
                + (imj=0)*(1+6*lambda);
    enddo;
  enddo;
?
? Fill in last 2 rows of W by hand:
?
  set ij = n*n-n-1;
  set w(ij) = 1+5*lambda;
  set ij = ij+1;
  set w(ij) = -2*lambda;
  set ij = n*n;
  set w(ij) = 1+lambda;
  mat tau = ser(w"y);
endproc HPtrend;
?
? ========================================================================
?
end;
noprint;
freq q;
smpl 1947:1 1986:3;
load rgnp;  ? US real gnp in $1982, from Perron
1056.5 1063.2 1067.1 1080.0 1086.8 1106.1 1116.3 1125.5
1112.4 1105.9 1114.3 1103.3 1148.2 1181.0 1225.3 1260.2
1286.6 1320.4 1349.8 1356.0 1369.2 1365.9 1378.2 1406.8
1431.4 1444.9 1438.2 1426.6 1406.8 1401.2 1418.0 1438.8
1469.6 1485.7 1505.5 1518.7 1515.7 1522.6 1523.7 1540.6
1553.3 1552.4 1561.5 1537.3 1506.1 1514.2 1550.0 1586.7
1606.4 1637.0 1629.5 1643.4 1671.6 1666.8 1668.4 1654.1
1671.3 1692.1 1716.3 1754.9 1777.9 1796.4 1813.1 1810.1
1834.6 1860.0 1892.5 1906.1 1948.7 1965.4 1985.2 1993.7
2036.9 2066.4 2099.3 2147.6 2190.1 2195.8 2218.3 2229.2
2241.8 2255.2 2287.7 2300.6 2327.3 2366.9 2385.3 2383.0
2416.5 2419.8 2433.2 2423.5 2408.6 2406.5 2435.8 2413.8
2478.6 2478.4 2491.1 2491.0 2545.6 2595.1 2622.1 2671.3
2734.0 2741.0 2738.3 2762.8 2747.4 2755.2 2719.3 2695.4
2642.7 2669.6 2714.9 2752.7 2804.4 2816.9 2828.6 2856.8
2896.0 2942.7 3001.8 2994.1 3020.5 3115.9 3142.6 3181.6
3181.7 3178.7 3207.4 3201.3 3233.4 3157.0 3159.1 3199.2
3261.1 3250.2 3264.6 3219.0 3170.4 3179.9 3154.5 3159.3
3186.6 3258.3 3306.4 3365.1 3444.7 3487.1 3507.4 3520.4
3547.0 3567.6 3603.8 3622.3 3655.9 3661.4 3683.3