' A simple RAS algorithm to estimate matrix cells when we know
' row and column totals

wfcreate a 2000 2006


' Matrix dimension
!md = 7
' Matrices used in the subroutine are global
MATRIX(!md,!md) Z1
matrix(!md,!md) z2
matrix(!md,!md) z3
matrix(!md) tr
matrix(!md) tc
matrix(!md,!md) r
matrix(!md,!md) s
matrix(!md) u
u.fill 1, 1, 1, 1, 1, 1
scalar sum

' Tolerance for conversion
scalar eps = 0.001

' Subroutine
subroutine ras
' requires global vector TC to hold column totals
' requires global vector TR to hold row totals

' Number of rows
!nr = !md
' Number of columns
!nc = !md

' Calculates matrix sum from column totals
sum = 0
for !i = 1 to !nc
   sum = sum+tc(!i)
next

' Initial values for Z1 matrix to be given by row totals
' assigned to cells as proportion of column totals over matrix total
for !i = 1 to !nr
   for !j = 1 to !nc

' Since we are working with consolidated figures we want the cells on the main
' diagonal to have zero values. Drop the IF statement if cells on the main diagonal
' can have values different than zero

      if (!i = !j) then
         z1(!i,!j) = 0
      else
        z1(!i,!j) = tr(!i)*tc(!j)/sum

      endif

   next
next

!converge=0
while !converge=0
   for !i = 1 to !nr
      sum=0
      for !j = 1 to !nc
       sum = sum+z1(!i, !j)
      next
     if (sum = 0) or (!i = !j) then
        r(!i,!i) = 0
     else
         r(!i,!i) = tr(!i)/sum
      endif
   next
   z2=r*z1
   
   for !j = 1 to !nc
      sum=0
      for !i = 1 to !nr
       sum = sum+z2(!i, !j)
      next
      if (sum = 0) or (!i = !j) then
         s(!j,!j) = 0
      else
         s(!j,!j) = tc(!j)/sum
      endif
   next
   z3 = z2*s
   !convcrit = 0
   for !i = 1 to !nr
      for !j = 1 to !nc
         if ((z3(!i,!j) - z1(!i,!j)) > eps) then
            !convcrit = !convcrit +1
         endif
      next
   next
   if !convcrit = 0 then
      !converge = 1
   else
      z1=z3
   endif
wend

endsub
' END OF SUBROUTINE
' Matrix Z3 holds results

' NOTE: if the procedure does not converge it will proceed indefinitely! In such a case, 
' use the Esc key to break out of the infinite loop


' As an example we use dividends paid and received by sector
' from the NIPA annual macroeconomic accounts
' for the years 2003 to 2006

' Dividends received are recorded in the rows
' Dividends paid go in columns

' starts reading data
series dh_r
dh_r.label(d) Households: Dividends received
dh_r.fill(o=2003) 422.6, 537, 598.9, 696.3
series dh_p = 0
dh_p.label(d) Households: Dividends paid
series dn_r = 0
dn_r.label(d) Non financial non corporate business: Dividends received
series dn_p = 0
dn_p.label(d) Non financial non corporate business: Dividends paid
series dc_r
dc_r.label(d) Non financial corporate business: Dividends received
dc_r.fill(o=2003) 55.8, 63.1, 307.2, 116.6
series dc_p
dc_p.label(d) Non financial corporate business: Dividends paid
dc_p.fill(o=2003) 348.5, 430.1, 506.4, 565.3
series df_r
df_r.label(d) Financial business: Dividends received
df_r.fill(o=2003) 77.3, 101.5, 134.8, 171.1
series df_p
df_p.label(d) Financial business: Dividends paid
df_p.fill(o=2003) 177.4, 226.3, 281.1, 345.5
series dg_r = 0
dg_r.label(d) Federal government: Dividends received
series dg_p = 0
dg_p.label(d) Federal government: Dividends paid
series ds_r
ds_r.label(d) State and Local government: Dividends received
ds_r.fill(o=2003) 2.2, 2.4, 2.4, 2.6
series ds_p = 0
ds_p.label(d) State and Local government: Dividends paid
series dw_r
dw_r.label(d) Rest of the World: Dividends received
dw_r.fill(o=2003) 68.9, 73.3, 92.5, 91.4
series dw_p
dw_p.label(d) Rest of the World: Dividends paid
dw_p.fill(o=2003) 100.9, 121.1, 348.3, 167.2

' Routine to fill the matrix

for %year 2003 2004 2005 2006
   ' Repeat for each year in the sample
   smpl %year %year
   ' Counter to fill the vectors with appropriate variables
   !count = 0
   ' Our sectors are identified by a single character:
   '   h > households
   '   n > Non financial non corporate
   '   c > Non financial corporate
   '   f > Financial
   '   g > Federal Government
   '   s > State and Local government
   '   w > Rest of the world
   for %s1 h n c f g s w
      !count = !count + 1
      tr(!count) = @elem(d{%s1}_r, %year)
      tc(!count) = @elem(d{%s1}_p, %year)
   next
   call ras
   !i = 0
   for %s1 h n c f g s w
      !i = !i + 1
      !j = 0
      for %s2 h n c f g s w
         !j = !j+1
' The resulting variables are labeled d_X_Y
' and record dividends paid from sector X to sector Y
         genr d_{%s2}_{%s1} = z3(!i,!j)
      next
   next
next

smpl @all