*
*  CONSUMER.PRG
*  Manual example 7.3
*
cal 1922
all 1962:1
open data consumer.wks
data(format=wks,org=cols) / wqfood wqvice wqdura wqrem dpfood dpvice dpdura dprem
set dq = wqfood + wqvice + wqdura + wqrem
nonlin(parmset=base) mu1 mu2 mu3 $
   p11 p12 p13 p21 p22 p23 p31 p32 p33
nonlin(parmset=relax) ep14 ep24 ep34
nonlin(parmset=symmetry) p12=p21 p13=p31 p23=p32
*
frml ffood wqfood = mu1*dq+p11*(dpfood-dprem)+$
   p12*(dpvice-dprem)+p13*(dpdura-dprem)+ep14*dprem
frml fvice wqvice = mu2*dq+p21*(dpfood-dprem)+$
   p22*(dpvice-dprem)+p23*(dpdura-dprem)+ep24*dprem
frml fdura wqdura = mu3*dq+p31*(dpfood-dprem)+$
   p32*(dpvice-dprem)+p33*(dpdura-dprem)+ep34*dprem
*
compute mu1=mu2=mu3=0.0
compute p11=p12=p13=p21=p22=p23=p31=p32=p33=0.0
compute ep14=ep24=ep34=0.0
*
nlsystem(parmset=base,sigma) / ffood fvice fdura
nlsystem(parmset=base+relax) / ffood fvice fdura
compute ep14=ep24=ep34=0.0
nlsystem(parmset=base+symmetry) / ffood fvice fdura
*
*  Additional code for non-linear maximum likelihood from
*  section 7.10.
*
dec frml[vect] ufrml
frml ufrml = ||wqfood-ffood,wqvice-fvice,wqdura-fdura||
dec symm sigma(3,3)
nonlin(parmset=sigmaparms) sigma
frml mvlikely = %logdensity(sigma,ufrml)
compute sigma=%identity(3)
maximize(parmset=base+sigmaparms,iters=100) mvlikely