@ ACH 9/11/06 @ 
@ Adaptive Conditional Heteroskedasticity, aka local Scale Model @
new;
library pgraph;
graphset;
_pdate = "";
_pltype = 6;
_pnum = 2;
_plwidth = 3;
_pgrid = {1,0};
_pmcolor=0;
_plwidth = 6;
_pcolor = {9, 10, 12};
zip = {};
missing = miss(1,1);
format /rdn 12,7;

load x[] = c:\data\crsp\vwd03.txt;
@ CRSP VWRET with Distributions, 12/25(first missing) - 12/2003 
  Arithmetic, from prior month.  Via Wharton Research Data Services. 
  936 obs. @
x = reshape(x,rows(x)/2,2);
ret = x[2:rows(x),2];  @ Delete missing 12/25 obs. @ 
n = rows(ret);
"n = "; n;
t = seqa(1926.0+1/12,1/12, n);
title ("CRSP value-weighted arithmetic returns with distributions");
_pline = 1~6~0~0~3000~0~1~0~6;
xy(t,ret);
"waiting"; wait;
ret = ln(1+ret);
load y[] = c:\data\crsp\Fama1moB.txt;
@ CRSP/FAMA 1 mo risk free rate, bid&ask (missing asks filled in),
  12/25 - 12/2004.  Annualized %, to subsequent month.  
  365 day year, continuously compounded. @ 
y = reshape(y, rows(y)/3,3);
rf = (y[.,2]+y[.,3])/2;
tt = seqa(1926.0, 1/12, rows(y));
title ("Fama 1 mo T-bill rate: bid, asked, avg.");
xy(tt,y[.,2:3]~rf);
"waiting"; wait;
rf = rf[1:n];
rf = rf/1200;
xsr = ret -rf;
_pcolor = {9, 12, 10};
title ("Log returns, lagged 1 mo Rf");
xy(t, ret~rf);
"waiting with ret, rf "; wait;

@ try arithmetic XS returns 
xsr = exp(xsr)-1; @

Title ("Excess log returns");
ytics(-.4,.4,.1,1);
xy(t, xsr);
"waiting"; wait;
{muml, semuml, vlogeta, h, v, d, dinf, u, w, lambda, gam, lr, 
   xbar, sexbar, muwls, sewls, mua, pdf, cdf} = ach(xsr);
format /rdn 12,7;
"muml = " muml;
"semuml " semuml;
"vlogeta" vlogeta; 
"d[n] = " d[rows(d)];
"dinf   " dinf;
"lr     " lr;
"xbar = " xbar;
"sexbar " sexbar;
"muwls  " muwls;
"sewls  " sewls;
" ";

for k(1,10,1);
ney = neyman(w[2:n],k);
"neyman " "k " ney cdfchic(ney,k);
endfor; 

_pcolor = {9,12};
_pcolor = {12,9};
title ("Demeaned excess returns, local scale");
_plegctl = {1 5 1955 .18};
_plegstr = "Demeaned XS returns\000Local scale";
xy(t, (xsr-muml)~sqrt(h));

"IGARCH estimates";
{mugarch, semu, omega, garchlam, garchgam, higarch, logl, lr} = 
   garch(xsr,1,.001,1);
"waiting"; wait;
format /rdn 12,7;
"mu     " mugarch;
"omega  " omega;
"lambda " garchlam;
"gamma  " garchgam;
"logl   " logl;
"lr     " lr;
title ("IGARCH demeaned excess returns, IGARCH local s.d.");
_plegstr = "IGARCH demeaned XS returns\000IGARCH local s.d.";
xy(t, (xsr-mugarch)~higarch);
"waiting"; wait;
title ("LSM local scale, IGARCH local s.d.");
ytics(0, .2, .05, 1);
_plegstr = "LSM local scale\000IGARCH local s.d.";
_plegctl = {1 5 1955 .12};
_pcolor = {9,12};
xy(t, sqrt(h)~higarch);
"waiting"; wait;

dcrop = d.*(d.>2) + 2.00000001.*(d.<=2);
lsmsd = sqrt(h.*dcrop./(dcrop-2));
title ("LSM local scale, LSM local s.d.");
_plegstr = "LSM local scale\000LSM local s.d.";
xy(t, sqrt(h)~lsmsd);
"waiting"; wait;
title ("LSM local s.d., IGARCH local s.d.");
_plegstr = "LSM local s.d.\000IGARCH local s.d.";
_pcolor = {9, 12};
xy(t, lsmsd~higarch);
"waiting"; wait;

"GARCH estimates";
{mugarch, semu, omega, garchlam, garchgam, hgarch, logl, lr} = 
   garch(xsr,0,.001,1);
"waiting"; wait;
format /rdn 12,7;
"mu     " mugarch;
"omega  " omega;
"lambda " garchlam;
"gamma  " garchgam;
"logl   " logl;
"lr     " lr;
title ("GARCH demeaned excess returns, GARCH local s.d.");
_plegstr = "GARCH demeaned XS returns\000GARCH local s.d.";
xy(t, (xsr-mugarch)~hgarch);
"waiting"; wait;
title ("LSM local scale, GARCH local s.d.");
ytics(0, .2, .05, 1);
_plegstr = "LSM local scale\000GARCH local s.d.";
_plegctl = {1 5 1955 .12};
_pcolor = {9,12, 10};
xy(t, sqrt(h)~hgarch);
"waiting"; wait;

title ("LSM local scale, IGARCH local s.d., GARCH local s.d.");
ytics(0, .2, .05, 1);
_plegstr = "LSM local scale\000IGARCH local s.d.\000GARCH local s.d.";
_plegctl = {1 5 1955 .12};
_pcolor = {9,12, 10};
xy(t, sqrt(h)~higarch~hgarch);
"waiting"; wait;

title ("Scale-adjusted errors");
ytics(-8,8,2,2);
_plegctl = 0;
_pcolor = {9, 12};
xy(t,u);
"waiting"; wait;
title ("Errors transformed to uniformity");
_pstype = 7;
_plctrl = -1;
_plegctl = 0;
_pcolor = 9;
ylabel ("Unit interval");
xy(t, w);
"waiting"; wait;
_plctrl = 0;
title ("Predictive DOF d(t)");
dshep = zeros(n,1);
kshep = zeros(n,1);
gshep = miss(1,1)*ones(n,1);
lamshep = miss(1,1)*ones(n,1);
omega = dinf/(dinf+1);
for i(2,n,1);
  dshep[i] = omega*(dshep[i-1]+1);
  kshep[i] = exp(digamma((dshep[i-1]+1)/2) - digamma(dshep[i]/2));
  gshep[i] = 1/(kshep[i]*dshep[i]);  
  lamshep[i] = dshep[i-1]./(kshep[i].*dshep[i]);
  endfor;
_plegctl = {1 5 11.5 1.5};
_plegstr = "d(t)\000d(inf)\000t-1\000Shephard d(t)";
_pltype = {6,1,5,3};
_pcolor = {9, 12, 12, 12};
xlabel("t");
ylabel ("");
ytics(0,10, 1, 1);
xtics(1,20,1,1);
xy(seqa(1,1,n), d~dinf*ones(n,1)~seqa(0,1,n)~dshep);
"waiting with d"; wait;

Title ("Gain gamma(t), attrition lambda(t), sum, Shephard gamma(t)");
d[1] = miss(1,1);
d[2] = miss(1,1);
d[3] = miss(1,1);
_plegctl = {1 5 11.5 0.3};
_plegstr = "gamma(t)\000lambda(t)\000sum\000Shephard gamma(t)";
_pltype = {6,6,6,3};
_pcolor = {9,10,12,9,12};
gg = gam.*d./(d-2);
ytics(0,1.2,.2,2);
xy(seqa(1,1,n),gam~lambda~(gam+lambda)~gshep);
"waiting"; wait;
xy(seqa(1,1,n),gam~lambda~(gam+lambda)~gshep~lamshep~(gshep+lamshep));
"waiting"; wait;
title ("");
ytics(0, 2, .2,2);
xy(seqa(1,1,n),gam~lambda~(gam+lambda)~gg~(gg+lambda));
"waiting"; wait;

title ("Histogram of transformed errors");
ylabel ("Frequency");
xlabel ("Unit interval");
ef = (n-1)/10;
_pline = 1~6~0~ef~3000~ef~1~0~6;
xtics(0,1,.1,1);
ytics(0,150,25,1);
xlabel ("Boundary");
{b,m,freq} = hist(w[2:n],seqa(.1,.1,10));
"waiting with hist "; wait;


title ("PP plot of transformed errors");
ylabel ("Empirical CDF"); 
xlabel ("Unit interval");
_plegctl = {1 5 .55 .15};
_plegstr = "Empirical CDF\000Theoretical CDF";
_pcolor = {9, 12};
xtics(0,1,.1,2);
ytics(0,1,.1,2);
ww = sortc(w[2:n],1);
pp = seqa(.5,1,n-1)/(n-1);
ww = 0|ww|1;
pp = 0|pp|1;
xy(ww,pp~ww);
end;

proc neyman(u,k);
local n, s, i, lm;
n = rows(u);
s = zeros(k,1);
i = zeros(k,k);
for j(1,k,1);
s[j] = sumc(u.^j) - n/(j+1);
for jj(1,k,1);
i[j,jj] = 1/(j+jj+1) - 1/((j+1)*(jj+1));
endfor;
endfor;
i = n*i;
lm = s'*inv(i)*s;
retp(lm);
endp;



proc(19) = ach(x);
@ Adaptive Conditional Heteroskedasticity, or Studentized IGARCH @       
@ x is nX1 time series with constant mean mu and stochastic precision @ 
@ returns {muml, semuml, vlogeta, h, v, d, dinf, u, w, lambda, lr, xbar, sexbar, muwls, sewls, mua, pdf, cdf} 
    muml = scalar ML est of mu
    vlogeta = var(log(eta)) = scalar variance of log of gamma-reverting precision shocks eta. 
    h = nX1 time-varying predictive harmonic variance = 1/E(precision of e[t]|e[t-1])
    d = nX1 DOF of predictive density for e[t]|e[t-1]
    dinf = d[inf]
    v = nX1 variance E(e[t]^2|e[t-1]) = h.*d./(d-2).  = inf for i = 1, 2, 3. 
    u = nX1 vector of e[t]/sqrt(h[t])  u[1] = missing code.
      u[t] should be T(u[t]; d[t]);
    w = nX1 vector of T^-1(u; d).  Should be U(0,1). w[1] = missing code.
    lambda = nX1 vector of GARCH(1,1)-like coefficients on v[t-1] for v[t]
    gam = nX1 vector of GARCH(1,1)-like coefficients on e[t-1]^2 for v[t] ("gamma")
    lr = Likelihood Ratio for v(log(eta)) = 0
    xbar, sexbar = OLS estimates
    muwls, sewls = WLS estimates using v
    mua = array of mu values for PDF, CDF
    PDF, CDF = posterior distribution of mu|x 
      (PDF inactivated, CDF not implemented) @
declare xglobach_; 
declare achdglob_, achhglob_, achuglob_, achlambdaglob_;
declare achgammaglob_;
local n, xbar, e, sexbar, p0, dp0, n1, p, loglik, its, loglik0, lr;
local muml, semuml,  vlogeta, h, v, d, u, w, muwls, sewls, mua, pdf, cdf;
local rz, iz, lambda, gam;
local s2wls, ewls, dmu, loglik1, loglik2, mumin, mumax, pdf2, simp, dinf;
xglobach_ = x;
n = rows(x);
xbar = sumc(x)/n;
e = x - xbar;
sexbar = sqrt(e'*e/(n*(n-1)));
"xbar = " xbar;
"se   = " sexbar; 
n1 = 1/n;
p0 = 0~xbar;
dp0 = n1~sexbar;
"Starting ACH search";
loglik0 = achloglik(p0);
@ veta bounded above by 10 to prevent possible anomalies @ 
{p, loglik, its} = rblob2(p0, dp0, &achloglik, .001, 1, 0, 10);
"Search ended";
muml = p[2];
vlogeta = p[1];
lr = 2*(loglik - loglik0);
h = achhglob_;
d = achdglob_;
u = achuglob_; 
w = 1-cdftc(u[2:n],d[2:n]);
w = miss(1,1)|w;
v = h.*d./(d-2);
lambda = achlambdaglob_;
gam = achgammaglob_;
rz = -1/0|2;
iz = indexcat(d,rz);
v[1] = 1/0;
v[iz] = v[1]*ones(rows(iz),1); 
sewls = sumc(1./v);
muwls = sumc(x./v)/sewls;
sewls = sqrt(1/sewls);
ewls = x-muwls;
s2wls = (ewls./v)'*ewls/(n-1);
"s2wls    " s2wls;
dmu = .01*sexbar;
p[2] = muml + dmu;
loglik2 = achloglik(p);
p[2] = muml - dmu;
loglik1 = achloglik(p);
semuml = (loglik2 + loglik1 - 2*loglik)/(dmu)^2;
semuml = sqrt(-1/semuml);
pdf = zeros(1001,1);
mua = zeros(1001,1);
cdf = 0;
/*
@ Don't bother with Posterior PDF for now @ 
"Computing Posterior PDF";
mumin = minc(muml|muwls)-5*semuml;
mumax = maxc(muml|muwls)+5*semuml;
dmu = (mumax-mumin)/1000;
mua = seqa(mumin, dmu, 1001);
for i(1,1001,1);
if i == 100*trunc(i/100); "i = " i; endif;
p[2] = mua[i];
pdf[i] = achloglik(p);
endfor;
"Done with Posterior PDF";
Title ("Log Posterior Density");
pdf = pdf - maxc(pdf);
xy(mua,pdf);
"waiting"; wait;
Title ("Posterior Density X Constant");
pdf = exp(pdf);
xy(mua,pdf);
"waiting"; wait;
Title ("Posterior Density");
pdf2 = reshape(pdf,501,2);
simp = dmu*(pdf[1] + pdf[1001] + 2*sumc(pdf2[2:500,1]) + 4*sumc(pdf2[1:500,2]))/3;
pdf = pdf/simp;
xy(mua, pdf);
"waiting"; wait;
*/
if vlogeta == 10;
  "Warning -- V(log(eta)) <= 10 is binding.";
  endif;
dinf = 2*apfromv(vlogeta,.0001,.0001);
retp(muml, semuml, vlogeta, h, v, d, dinf, u, w, lambda, gam, lr, xbar, sexbar, muwls, sewls, mua, pdf, cdf);
endp;

proc achloglik(p);
@ p = vlogeta~mu @ 
@ receives x as global xglobach_ @
@ vlogeta >= 0 is variance of log of eta-shocks to gamma distributed 
  precision posterior.  @ 
@ returns loglik directly, d and v as globals achdglob_ and achvglob_ @ 
local e, ee, loglik, s, h, vlogeta, mu, u, ap, bp, lambda, n, k;
local apmax,gam;
vlogeta = p[1];
mu = p[2];
e = xglobach_ - mu;
ee = e.*e;
n = rows(e);
ap = zeros(n,1);
bp = zeros(n,1);
u = zeros(n,1);
h = zeros(n,1);
lambda = zeros(n,1);
gam = zeros(n,1);
u[1] = miss(1,1);
h[1] = miss(1,1);
lambda[1] = miss(1,1);
gam[1] = miss(1,1);
loglik = 0;
if vlogeta < 2/n^2;
 apmax = n/2;
 else; 
 apmax = apfromv(vlogeta, .0001, .0001);
 endif;
for i(2,n,1);  
ap[i] = trigaminv(trigamma(ap[i-1]+.5)+vlogeta, apmax, .001, .001);
k = exp(digamma(ap[i-1] + .5) - digamma(ap[i]));
bp[i] = (bp[i-1] + ee[i-1]/2)/k;
h[i] = bp[i] / ap[i];
s = sqrt(h[i]);
lambda[i] = ap[i-1]/(k*ap[i]);
gam[i] = 1/(k*2*ap[i]);
u[i] = e[i]/s;
loglik = loglik + lnpdft(u[i], 2*ap[i]) - ln(s);
endfor;
achdglob_ = 2*ap;
achhglob_ = h;
achuglob_ = u;
achlambdaglob_ = lambda;
achgammaglob_ = gam;
retp(loglik);
endp;

proc v0loglik(mu);
retp(achloglik(0~mu));
endp;

proc lnpdft(x,d);
@ x is n X 1 vector @
@ d is n X 1 vector of DOFs @ 
@ returns n X 1  t(d) log PDF @
@ note that each x[i] may have a different d[i] @ 
retp(lngam((d+1)/2) - lngam(d/2) - ln(d*pi)/2 - ((d+1)/2).*ln(1+(x.^2)./d) );
endp;

proc lngam(z);
@ Computes ln(gamma(z)) for real positive z , using 
  lnfact for z < 7 and ln(gamma) for z <= 7.  
  ln(gamma(z)) = lnfact(z-1).  Howver, lnfact is only defined 
  for positive arguments, and does worse for z < 7.  
  gamma(z) actually stops growing for z > 160 or so @   
local zz, zzz;
zz = z + (z .< 1);    @ zz = z for z >= 7 @
zzz = minc(z'|10*ones(1,rows(z)));  @ zzz = z for z < 7 @
retp(lnfact(zz-1).*(z.>=7) + ln(gamma(zzz)).*(z.<7)); 
endp;

proc apfromv(vlogeta,da,dv);
@ returns aprime, within da, such that 
    vlogeta = trigamma(aprime) - trigamma(aprime + .5), 
    to within dv.
  vlogeta is a scalar.  @ 
local x1, y1, x2, y2, xx, yy, y;
y = vlogeta;
x1 = 1/sqrt(2*vlogeta);
y1 = trigamma(x1) - trigamma(x1+.5);
if y1 > y;
 double:
 x2 = 2*x1;
 y2 = trigamma(x2) - trigamma(x2+.5);
 if y2 < y; 
   goto split;
 else;
   x1 = x2;
   y1 = y2;
   goto double;
 endif;
else;
  halve:
  x2 = x1;
  y2 = y1;
  x1 = x2/2;
  y1 = trigamma(x1) - trigamma(x1+.5);
  if y1 > y;
    goto split;
  else;
    goto halve;
  endif;
endif;
split: 
if (x2-x1)<da and (y1-y2)<dv;
 retp((x1+x2)/2);
endif;
xx = (x1+x2)/2;
yy = trigamma(xx) - trigamma(xx+.5);
if yy > y;
 x1 = xx;
 y1 = yy;
 goto split;
else;
 x2 = xx;
 y2 = yy;
 goto split;
endif;
endp;


proc trigaminv(y,x0,dx,dy);
@ finds x, within dx, such that trigamma(x) = y, to within dy.
  x0 is initial guess.  
  y is a scalar.  @ 
local x1, y1, x2, y2, xx, yy;
x1 = x0;
y1 = trigamma(x1);
if y1 > y;
 double:
 x2 = 2*x1;
 y2 = trigamma(x2);
 if y2 < y; 
   goto split;
 else;
   x1 = x2;
   y1 = y2;
   goto double;
 endif;
else;
  halve:
  x2 = x1;
  y2 = y1;
  x1 = x2/2;
  y1 = trigamma(x1);
  if y1 > y;
    goto split;
  else;
    goto halve;
  endif;
endif;
split: 
if (x2-x1)<dx and (y1-y2)<dy;
 retp((x1+x2)/2);
endif;
xx = (x1+x2)/2;
yy = trigamma(xx);
if yy > y;
 x1 = xx;
 y1 = yy;
 goto split;
else;
 x2 = xx;
 y2 = yy;
 goto split;
endif;
endp;