*
* @ChowLin   oldser  start  end  newser
* # list of related series in regression format.
*
* ChowLin distributes a series, changing the frequency to a higher one while maintaining
* the sum over each period. The newer DISAGGREGATE.SRC procedure does the same calculation,
* and more.
*
* Parameters:
*
*   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.
*
*   start end  = range over which to do the calculations. By default, the defined
*            range of <<oldser>>
*
*   newser = output 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  =input value of the AR1 parameter for AR1 and RWAR1 models
*   ESTIMATE/[NOESTIMATE]
*     If ESTIMATE, the RHO is estimated (by maximum likelihood) for the AR1 and
*     RWAR1 models.
*
*   FACTOR=Increase in the recording frequency (3 for quarterly to monthly)
*
*
* Notes:
*   Make sure that you include the CONSTANT if you're using MODEL=AR1. You <<cannot>> include
*   the CONSTANT for any of the other models.
*
* Revision Schedule:
*   01/2005 Written by Tom Doan, Estima
*   04/2005 Normalize startl, test estimation of rho
*   04/2007 Correct multiple declaration of FACTOR option
*   03/2008 Revisions to comments
*
* References:
*
*   MODEL=AR1 is Chow and Lin(1971), "Best Linear Unbiased Interpolation,
*     Distribution and Extrapolation of Time Series by Related Series", Review of
*     Economics and Statistics, vol 53, pp 372-375.
*
*   MODEL=RW is Fernandez(1981), "A Methodological Note on the Estimation of Time
*     Series", Review of Economics and Statistics, vol 63, pp 471-478.
*
*   MODEL=RWAR1 is Litterman(1983), "A Random Walk, Markov Model for the
*     Distribution of Time Series", JBES, vol 1, pp 169-173, with a different handling
*     of the initial periods.
*
* Variables Defined:
*   %LOGL   log likelihood
*
procedure ChowLin y start end ydist
type series  y *ydist
type integer start end
*
option integer factor 4
option real    rho    .9
option choice  model  1  rw1  ar1  rwar1  rw2
option switch  estimate
*
local integer i j n
local integer startl endl
local equation xvar
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,print,exact,g=g,var=concentrate,$
      a=af,c=cf,y=yf,sw=sw,reject=abs(rhox)>=1.00) startl endl xstates
}
else
   dlm(type=smooth,free=n,method=solve,print,exact,g=g,var=concentrate,$
     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


