set echo off set messages off open hpfilter.gdt --quiet scalar DEBUG = 0 scalar lambda = 6.25 smpl 1962 2010 scalar n=\$nobs series yseries=DLGDP series tauhat sess seas matrix y = {DLGDP} # Construct the K matrix matrix K = zeros(n-2,n) loop i=1..(n-2) --quiet K[i,i:i+2] = {1, -2, 1} end loop matrix A = I(n) + K'K .* lambda matrix invA = invpd(A) series tauhat = invA*y series cycle=yseries-tauhat if DEBUG printf "check: %g\n", sum(abs(cycle - hpfilt(yseries, lambda))) endif /* for simplicity, first assume that V(y) = scalar matrix */ matrix var_tauhats = qform(invA, var(yseries) .* I(n)) series sess = sqrt(diag(var_tauhats)) series lowers=tauhat-1.96*sess series uppers=tauhat+1.96*sess gnuplot tauhat lowers uppers --time-series --output=display --with-lines /* now let's be more realistic and specify V(y) more appropriately. I've identified an AR(1) process for the series. */ # construct V(y) if DEBUG ols yseries const yseries(-1) endif scalar ssq=0.000392 scalar rho=0.411384 matrix va = rho .^ abs(seq(1,n) .- seq(1,n)') matrix va = va .* (ssq/(1-rho^2)) series seas = sqrt(diag(qform(invA,va))) series lowera=tauhat-1.96*seas series uppera=tauhat+1.96*seas gnuplot tauhat lowera lowers uppera uppers --time-series --output=display --with-lines