*
* ChowLin distributes a series, changing the frequency to
* a higher one while maintaining the sum over each period
*
* Syntax:
*
*  @ChowLin   oldser  start  end  newser
*  # list of related series in regression format. (Make sure that
*    you include the CONSTANT).
*
*    oldser = series to distribute. This should be set up in the higher
*             frequency.  ChowLin uses the final subperiod as the value
*             to be distributed.
*
*             If you read the series from a RATS format file in the higher
*             frequency, it will be set up this way automatically.
*
*    newser = higher frequency series. Note that, if you want the output
*             to be at annual rates, you'll have to scale this up at the
*             end, because ChowLin distributes the sum, not the average.
*
*  Options:
*     MODEL=[RW1]/AR1/RWAR1/RW2
*       Specifies one of several statistical models for the error process.
*       RW1 is a random walk, AR1 is a first order autoregression,
*       RWAR1 is an ARIMA(1,1,0) and RW2 is ARIMA(0,2,0).
*
*     RHO  =value of the AR1 parameter for AR1 and RWAR1 models
*
*     FACTOR=Increase in the recording frequency (3 for quarterly to monthly)
*
*  Revision Schedule
*   01/2005 Written by Tom Doan, Estima
*   04/2005 Normalize startl, test estimation of rho
*
procedure ChowLin y start end ydist
type series  y *ydist
type integer start end
option integer factor 4
option real    rho    .9
option integer factor 3
option choice  model  1  rw1  ar1  rwar1  rw2
option switch  estimate
*
local integer i j n
local integer startl endl
local equation xvar
local rect a c
local symm sw
local frml[rect] af cf
local frml yf
local equation xvar
local rect abase
local series[vect] xsums xstates
local vect ones sumv
local rect g
local real rhox
*
* Get the range, taking the range of the input series as the default
* Adjust the start date to the beginning of the first period
*
inquire(series=y) startl<<start endl<<end
if %clock(startl,factor)<>1
   compute startl=startl-%clock(startl,factor)+1
*
* Pull in the regressors, and generate a SERIES[VECT] of their sums over
* "factor" periods
*
equation xvar *
compute n=%eqnsize(xvar)
gset xsums startl endl = sumv=%zeros(n,1),%do(j,0,factor-1,sumv=sumv+%eqnxvector(xvar,t-j)),sumv
*
* The state vector will be the states required for the error process, augmented
* to cover lags up to factor-1, concatenated with the regression coefficients.
* The latter will be estimated using the diffuse Kalman smoother.
*
dim ones(factor)
compute ones=%const(1.0)
frml cf = %blockglue(||ones|xsums||)
*
* The only error term hitting the state vector will be for the first
* state equation.
*
compute sw = %outerxx(%unitv(n+factor,1))
*
* The "A" matrix will have the same form for all models except for the first
* row at the top. The next factor-1 rows will have the typical structure
* for lagged states, while the bottom corner will be the identity. Because
* the error states will have different stationarity properties for different
* models, the "G" matrix will also have different structures.
*
dim abase(n+factor,n+factor)
ewise abase(i,j)=%if(i<=factor,(i==j+1),(i==j))
if model==1 {
   frml af = abase(1,1)=1.0,abase
   dim g(factor-1,n+factor)
   fmatrix(diffs=1) g
}
else
if model==2 {
   frml af = abase(1,1)=rhox,abase
   dim g(factor,n+factor)
   ewise g(i,j)=i==j
}
else
if model==3 {
   frml af = abase(1,1)=1+rhox,abase(1,2)=-rhox,abase
   dim g(factor-1,n+factor)
   fmatrix(diffs=1) g
}
else
if model==4 {
   frml af = abase(1,1)=2,abase(1,2)=-1,abase
   dim g(factor-2,n+factor)
   fmatrix(diffs=2) g
}
*
* Only allow y to be observed in the final subperiod
*
frml yf = %if(%clock(t,factor)==factor,y,%na)
*
* Apply Kalman smoothing. The regression coefficients will be at
* the end, and the smoothed error terms at the higher frequency will
* be in the first component.
*
compute rhox=rho
if estimate {
   nonlin rhox
   dlm(type=smooth,free=n,method=bfgs,noprint,exact,g=g,scale,a=af,c=cf,y=yf,sw=sw) startl endl xstates
}
else
   dlm(type=smooth,free=n,method=solve,noprint,exact,g=g,scale,a=af,c=cf,y=yf,sw=sw) startl endl xstates

report(action=define)
report(atrow=1,atcol=1,span) 'Chow-Lin Distribution'
report(atrow=2,atcol=1,span) 'Series '+%l(y)
report(atrow=4,atcol=1) 'Label' 'Coefficient'
do i=1,n
   report(row=new,atcol=1) %eqnreglabels(xvar)(i) xstates(endl)(factor+i)
end do i
report(action=format,atcol=2,width=8)
report(action=show)
*
* Generate the distributed series
*
set ydist startl endl = xstates(t)(1)+%dot(%eqnxvector(xvar,t),%xsubvec(xstates(t),factor+1,factor+n))
end

