*
* @Disaggregate( options )  oldser  start  end  newser
* # list of related series in regression format.
*
* Does a temporal disaggregation of a series, changing the frequency to
* a higher one while maintaining the sum, average or final value over
* each period of the levels.
*
* Make sure that you include the CONSTANT if you're using TSMODEL=AR1. You
* <<cannot>> include the CONSTANT for any of the other models.
*
* Parameters:
*
* oldser = series to distribute. This should be set up in the higher
*          frequency. Distribute 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.
*
* Options:
*     MODEL=LINEAR/LOGLIN/MULTIPLICATIVE
*        Specifies how the time series model and the regression model combine
*        to model the observed data.
*         MODEL=LINEAR --> y=regression + noise.
*         MODEL=LOGLIN --> y=exp(regression+noise)
*         MODEL=MULT   --> y=regression x noise
*
*     MAINTAIN=[SUM]/AVERAGE/FINAL/WEIGHTEDSUM
*        Indicates what the disaggregated data are supposed to match in the
*        lower frequency data. SUM is the sum across the subperiods, AVERAGE is
*        the average, FINAL is the last value, and WEIGHTEDSUM is the weighted
*        sum with weights given by the WEIGHTS option.
*
*     TSMODEL=[RW1]/AR1/RWAR1/RW2
*       Specifies one of several statistical models for the error process.
*       RW1 is a random walk, AR1 is a (stationary) first order autoregression,
*       RWAR1 is an ARIMA(1,1,0) and RW2 is ARIMA(0,2,0).
*
*     WEIGHTS=series of weights for MAINTAIN=WEIGHTEDSUM
*
*     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)
*
*     PRINT/[NOPRINT]
*     TITLE=title for report ["Disaggregation Analysis"]
*
* References:
*
* MODEL=LINEAR,TSMODEL=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=LINEAR,TSMODEL=RW1 is Fernandez(1981), "A Methodological Note on the
*    Estimation of Time Series", Review of Economics and Statistics, vol 63, pp
*    471-478.
*
* MODEL=LINEAR,TSMODEL=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.
*
* MODEL=MULTIPLICATIVE,TSMODEL=RW1 is the proportional Denton method
*
* MODEL=LINEAR,TSMODEL=RW1,MAINTAIN=SUM with no related series is a
*    generalization of Boot, Feibes, Lisman(1967) "Further methods of
*    derivation of quarterly figures from annual data", Applied
*    Statistics, vol. 16, no 1, 65-75.
*
* Variables Defined:
*  %LOGL   log likelihood
*
* Revision Schedule:
*  04/2008 Written by Tom Doan, Estima
*  11/2010 Correct output labeling of the AR(1) and RW models.
*  03/2013 Add MAINTAIN=WEIGHTEDSUM. Correct several problems with
*          output.
*  07/2013 Allow for FACTOR=1 to patch missing data without changing
*          frequency.
*  08/2016 Correct MODEL=LOGLIN, MAINTAIN=FINAL
*
procedure Disaggregate y start end ydist
type series  y *ydist
type integer start end
*
option integer factor   4
option real    rho      .9
option choice  tsmodel  1  rw1    ar1     rwar1  rw2
option choice  model    1  linear loglin  multiplicative
option choice  maintain 1  sum    average final  weightedsum
option switch  estimate 1
option switch  print    0
option series  weights
option string  title
*
local integer      i j n nfree
local integer      iters itlimit
local integer      startl endl
local equation     xvar
local symm         sw
local frml[rect]   af
local frml[vect]   cf
local frml[vect]   lcf
local frml         yf ydistf
local rect         abase
local series[vect] xsums xstates xstate0
local series[symm] vstates
local series       observed
local vect         sumv mx xv xt cv beta lc
local rect         g
local real         rhox yadj
local string       modeldesc tsmodeldesc ltitle
local report       dreport
local integer      nrequired ndlm
*
* 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)
*
* The multiplicative model needs to peg one of the regression
* coefficients (to 1.0), so the number of free regression coefficients
* is the number of input variables - 1.
*
if model==3
   compute nfree=n-1
else
   compute nfree=n
dim xt(n)
*
* 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.
*
* The only error term hitting the state vector will be for the first
* state equation.
*
* This makes sure that we have at least the minimal required size for
* the state-space model of the error term.
*
if tsmodel==1.or.tsmodel==2
   compute nrequired=1
else
   compute nrequired=2
*
compute ndlm=%imax(factor,nrequired)

compute sw = %outerxx(%unitv(ndlm+nfree,1))
*
* The "A" matrix will have the same form for all models except for the
* first row at the top. The next ndlm-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(ndlm+nfree,ndlm+nfree)
ewise abase(i,j)=%if(i<=ndlm,(i==j+1),(i==j))
if tsmodel==1 {
   frml af = abase(1,1)=1.0,abase
   dim g(ndlm-1,ndlm+nfree)
   fmatrix(diffs=1) g
   compute tsmodeldesc="RW"
}
else
if tsmodel==2 {
   frml af = abase(1,1)=rhox,abase
   dim g(ndlm,ndlm+nfree)
   ewise g(i,j)=(i==j)
   compute tsmodeldesc="AR(1)"
}
else
if tsmodel==3 {
   frml af = abase(1,1)=1+rhox,abase(1,2)=-rhox,abase
   dim g(ndlm-1,ndlm+nfree)
   fmatrix(diffs=1) g
   compute tsmodeldesc="ARIMA(1,1,0)"
}
else
if tsmodel==4 {
   frml af = abase(1,1)=2,abase(1,2)=-1,abase
   dim g(ndlm-2,ndlm+nfree)
   fmatrix(diffs=2) g
   compute tsmodeldesc="ARIMA(0,2,0)"
}
*
* Only allow y to be observed in the final subperiod
*
set observed startl endl = %clock(%period(t),factor)==factor
*
* Set the linear combination of disaggregated values which map to the
* observed y
*
frml lcf = lc
dim lc(ndlm)
if maintain==1
   ewise lc(i)=%if(i<=factor,1.0,0.0)
else
if maintain==2
   ewise lc(i)=%if(i<=factor,1.0/factor,0.0)
else
if maintain==3
   ewise lc(i)=%if(i==1,1.0,0.0)
else {
   compute lc=%zeros(ndlm,1)
   if .not.%defined(weights)
      disp "### Need WEIGHT option for MAINTAIN=WEIGHTEDSUM"
   frml lcf = %do(i,1,factor,lc(i)=weights(t+1-i)),lc
}
*
dim cv(ndlm+nfree) mx(ndlm)
*
if model==1 {
   *
   * Linear model. No iterations are required, no initial expansion
   * point is necessary.
   *
   compute itlimit=1,modeldesc="Linear"
   if n>0 {
      gset xsums startl endl = lc=lcf,sumv=%zeros(n,1),%do(i,1,factor,sumv=sumv+lc(i)*%eqnxvector(xvar,t+1-i)),sumv
      frml cf = lcf(t)~~xsums
      frml ydistf = beta=%xsubvec(xstate0,ndlm+1,ndlm+n),xv=%eqnxvector(xvar,t),$
         xstate0(t)(1)+%dot(beta,xv)
   }
   else {
      frml cf = lcf
      frml ydistf = xstate0(t)(1)
   }
   frml yf = %if(%clock(t,factor)==factor,y,%na)
}
else
if model==2 {
   *
   * Log linear model. Start with an initial expansion point of
   * log(y*lc) for the first <<factor>> lags of the disturbance (which
   * will almost always cover all of those), and zeros for any other
   * required lags of the disturbances and for the beta coefficients.
   *
   compute itlimit=10,modeldesc="Log-Linear"
   *
   * 08/2016 change to not take log zero on MAINTAIN=FINAL
   *
   gset xstate0 startl endl = xv=%zeros(n+ndlm,1),lc=lcf,%do(i,1,factor,xv(i)=%if(lc(i)==0,0.0,log(y*lc(i)))),xv
   *
   * This generates the <<c>> formula for the linearized measurement equation.
   *
   if n>0 {
      frml cf = lc=lcf,beta=%xsubvec(xstate0,ndlm+1,ndlm+n),cv=%zeros(n+ndlm,1),$
       %do(i,1,factor,%(xv=%eqnxvector(xvar,t+1-i),cv(i)=mx(i)=lc(i)*exp(xstate0(t)(i)+%dot(beta,xv)),$
        %do(j,1,n,cv(ndlm+j)=cv(ndlm+j)+mx(i)*xv(j)))),$
         cv
      frml ydistf = beta=%xsubvec(xstate0,ndlm+1,ndlm+n),xv=%eqnxvector(xvar,t),$
       exp(xstate0(t)(1)+%dot(beta,xv))
   }
   else {
      frml cf     = lc=lcf,cv=%zeros(ndlm,1),%do(i,1,factor,%(cv(i)=mx(i)=lc(i)*exp(xstate0(t)(i)))),cv
      frml ydistf = exp(xstate0(t)(1))
   }
   *
   * This generates the <<y>> formula for the linearized measurement
   * equation. This is only observable in periods which are multiples of
   * <<factor>>
   *
   frml yf = %if(.not.observed,%na,$
      %(yadj=y+%dot(cf(t),xstate0(t)),%do(i,1,factor,yadj=yadj-mx(i)),yadj))
}
else
if model==3 {
   *
   * Multiplicative model. Start with an initial expansion point of lc *
   * y(t)/x(t)(1). The coefficient on the first X variable is pegged to
   * one in this, so we start with all other coefficients at zero.
   *
   * If there is only one input variable, there's no need to iterate
   *
   if nfree==0
      compute itlimit=1
   else
      compute itlimit=10
   compute modeldesc="Multiplicative"
   gset xstate0 startl endl = lc=lcf,xv=%zeros(ndlm+nfree,1),xt=%eqnxvector(xvar,t),$
      %do(i,1,factor,xv(i)=lc(i)*y/xt(1)),xv
   *
   * This generates the <<c>> formula for the linearized measurement equation.
   *
   if n==1 {
      frml cf = lc=lcf,cv=%zeros(ndlm,1),$
       %do(i,1,factor,%(xv=%eqnxvector(xvar,t+1-i),cv(i)=lc(i)*xv(1),mx(i)=cv(i)*xstate0(t)(i))),cv
      frml ydistf = xv=%eqnxvector(xvar,t),xstate0(t)(1)*xv(1)
   }
   else {
      frml cf = lc=lcf,beta=1.0~~%xsubvec(xstate0,ndlm+1,ndlm+nfree),cv=%zeros(ndlm+nfree,1),$
       %do(i,1,factor,%(xv=%eqnxvector(xvar,t+1-i),cv(i)=lc(i)*%dot(beta,xv),mx(i)=cv(i)*xstate0(t)(i),$
        %do(j,1,nfree,cv(ndlm+j)=cv(ndlm+j)+lc(i)*xv(j+1)*xstate0(t)(i)) ) ),$
         cv
      frml ydistf = beta=1.0~~%xsubvec(xstate0,ndlm+1,ndlm+nfree),xv=%eqnxvector(xvar,t),$
         xstate0(t)(1)*%dot(beta,xv)
   }
   *
   * This generates the <<y>> formula for the linearized measurement
   * equation. This is only observable in periods which are multiples of
   * <<factor>>
   *
   frml yf = %if(.not.observed,%na,$
      %(yadj=y+%dot(cf(t),xstate0(t)),%do(i,1,factor,yadj=yadj-mx(i)),yadj))
}
*
* 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
do iters=1,itlimit
   if estimate {
      nonlin rhox
      dlm(type=smooth,free=nfree,method=bfgs,noprint,exact,g=g,var=concentrate,$
         a=af,c=cf,y=yf,sw=sw,reject=abs(rhox)>=1.00) startl endl xstates vstates
   }
   else
      dlm(type=smooth,free=nfree,method=solve,noprint,exact,g=g,var=concentrate,$
        a=af,c=cf,y=yf,sw=sw) startl endl xstates vstates
   *
   * Reset the expansion point
   *
   gset xstate0 startl endl = xstates(t)
end do iters
if %defined(title)
   compute ltitle=title
else
   compute ltitle="Disaggregation Analysis"
report(use=dreport,action=define,title=ltitle)
report(use=dreport,atrow=1,atcol=1,span) ltitle
report(use=dreport,atrow=2,atcol=1,span) "Series: "+%l(y)
report(use=dreport,atrow=3,atcol=1,span) modeldesc+" with "+tsmodeldesc+" time series errors"
report(use=dreport,atrow=5,atcol=1) "Label" "Coefficient" "StdError"
compute %nreg=n
dim %beta(%nreg) %stderrs(%nreg) beta(nfree) xv(nfree)
ewise beta(i)=xstates(endl)(ndlm+i)
ewise xv(i)=sqrt(%variance)*vstates(endl)(ndlm+i,ndlm+i)
if nfree==n
   compute %beta=beta,%stderrs=xv
else
   compute %beta=1.0~~beta,%stderrs=0.0~~xv
do i=1,n
   report(use=dreport,row=new,atcol=1) %eqnreglabels(xvar)(i) %beta(i) %stderrs(i)
end do i
report(use=dreport,action=format,atcol=2,width=8)
if print
   report(use=dreport,action=show)
*
* Generate the distributed series
*
if %defined(ydist)
   set ydist startl endl = ydistf
end
