*
* @DISTRIB(options)   oldser   newser
*
* DISTRIB computes a distribution of a series, changing the frequency to
* a higher one while maintaining the sum each original period, e.g.
* producing a monthly "GNP" estimate from quarterly GNP. The newer
* procedure DISAGGREGATE.SRC can do the same calculation and quite a bit
* more.
*
* If your original series is at an annual rate or is indexed, and you
* wish to maintain that, the output series must be multiplied by the
* ratio between periodicities, e.g. 3 for quarterly to monthly.
*
* Parameters:
*   oldser = series to distribute.  This should be set up in the higher
*            frequency (monthly in the example).  DISTRIB takes as the value
*            to distribute that at the last month of each quarter.
*
*            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
*
* Options:
*   MODEL=[RW1]/AR1/RWAR1/RW2
*     Specifies one of several statistical models for the 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)
*
*
* Distribution Procedure, solved using DP algorithm
*
*    X(t)=(x(3,t),x(2,t),x(1,t),1)'
*    U(t)=(u(3,t),u(2,t),u(1,t))'
*
*    X(t)=A X(t-1) + B U(t)
*
*    Problem is to minimize the sum of U(t)' R U(t) (R = identity),
*    subject to c(t)X(t) = 0, c(t)=(1,1,1,-y(t))
*
*    With cost-to-go functional  X(t)'Q(t)X(t), the solution is
*    of the form U(t)=D(t)X(t-1)
*
*    Total "Cost" is displayed and also is available as the variable
*    %FUNCVAL
*
* Revision Schedule
*    03/2005 Rewritten for version 6 by Tom Doan, Estima
*
* Reference:
*    The special case of MODEL=RW1 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.
*
procedure  distrib  oldser  newser
type series oldser
type series *newser
*
option  real     rho     .9
option  integer  factor  3
option  choice   model   1  rw1  ar1  rwar1  rw2
*
local rect a
local rect b
local rect c
local symm r
local rect v
local symm q
local rect d
local series[rect] dmats
local rect lambda lc
local symm temp1 qq
local vect qb
local vect y yshort
local real predict error scalar
local real rhol
*
local integer n1 n2 total nper i j
*
inquire(series=oldser) n1 n2
compute total=n2-n1+1
compute nper =total/factor
set newser n1 n2 = 0.0
gset dmats 1 nper = %zeros(factor,factor+1)
*
dim a(factor+1,factor+1) b(factor+1,factor)  c(1,factor+1)  r(factor,factor) $
    q(factor+1,factor+1)
dim qb(factor+1)
*
compute rhol=rho
if model==1.or.model==2
   {
   if model==1
      compute rhol=1.0
   dim v(factor,factor-1)
   ewise qb(i)=rhol**(i-1)
   ewise a(i,j)=(i <= factor .and. j==1) * qb(factor+2-i)
   ewise v(i,j)=(i==j)-rhol*(i==j+1)
   }
else
   {
   if model==4
      compute rhol=1.0
   dim v(factor,factor-2)
   compute qb(1)=1.0
   do i=2,factor+1
      compute qb(i)=rhol*qb(i-1)+1.0
   end do i
   compute a=%const(0.0)
   do i=1,factor
      compute a(i,1)=qb(factor+2-i)
      compute a(i,2)=-rhol*qb(factor+1-i)
   end do i
   ewise v(i,j)=(i==j)-(1+rhol)*(i==j+1)+rhol*(i==j+2)
   }
*
compute a(factor+1,factor+1)=1.0
ewise b(i,j)=%if( i <= j.and.i <= factor, qb(%iabs(i-j)+1) , 0.0)
compute c=%const(1.0)
compute r=%mscalar(1.0)
*
*  DP transition matrices are set up.  Loop over all periods except the
*  first.
*
compute q=%const(0.0)
*
do i=nper,2,-1
   compute c(1,factor+1)=-oldser(n1+i*factor-1)

   compute temp1=inv(r+%mqform(q,b))
   compute scalar=%qform(temp1,c*b)
   compute lambda=(1./scalar)*(c*b*temp1*tr(b)*q-c)
   compute d=temp1*(tr(b)*tr(c)*lambda-tr(b)*q)*a
   compute q=%mqform(q,a+b*d)+tr(d)*d
   compute dmats(i)=d
end do i
*
*  Compute minimizer for the initial period
*
overlay q(1,1)        with qq(factor,factor)
overlay q(factor+1,1) with qb(factor)
*
compute qq=qq+v*tr(v)
*
dim y(factor+1)
overlay c(1,1) with lc(1,factor)
*
compute predict=%scalar(lc*inv(qq)*qb)
compute error  =oldser(n1+factor-1)+predict
compute yshort =error*inv(qq)*tr(lc)*inv(lc*inv(qq)*tr(lc))-inv(qq)*qb
do i=1,factor
  compute newser(n1+factor-i)=y(i)=yshort(i)
end do j
compute y(factor+1)=1.0
compute %funcval=%qform(q,y)
display "Distribution criterion function = " %funcval
*
do i=2,nper
  compute d=dmats(i)
  compute y=(a+b*d)*y
  do j=1,factor
     compute newser(n1+i*factor-j)=y(j)
  end do j
end do i
end

