/*
   Interpolate quarterly values for an annual time-series,
   via Chow-Lin method. See Gregory C. Chow and An-loh Lin,
   "Best Linear Unbiased Interpolation, Distribution, and 
   Extrapolation of Time Series by Related Series",
   The Review of Economics and Statistics, Vol. 53, No. 4 
   (Nov., 1971) pp. 372-375.
*/

function scalar afunc (scalar a)
  scalar num = a^5 + 2*a^4 + 3*a^3 + 2*a^2 + a
  scalar den = 3 + 2*a^2 + 4*a
  return num/den
end function

function scalar chow_lin (matrix *ma, scalar targ)
  scalar val = afunc(ma[1])
  return -(val - targ)^2
end function

function matrix autocorr_matrix (scalar a, scalar n)
  matrix Vcol = zeros(n,1)
  scalar apow = 1
  scalar p = n
  scalar r = n*(n+1)/2
  
  # use vech approach to avoid n^2 genr calls

  loop i=1..n -q
    Vcol[i] = apow
    apow *= a
  endloop

  matrix vechV = zeros(r,1)
  r = 1
  p = n
  loop i=1..n -q
    vechV[r:r+p-1,] = Vcol[1:p,]
    r += p
    p--
  endloop
  
  return unvech(vechV)
end function

function matrix interpol_aq (series y_annual, scalar T)
  # number of quarterly observations
  scalar T4 = T * 4

  # constant and quadratic trend
  matrix Xq = ones(T4, 3)
  Xq[,2] = seq(1,T4)'
  Xq[,3] = Xq[,2].^2
  matrix Crow = zeros(T4+4,1)
  Crow[1:4,] = 1
  matrix C = mshape(Crow, T4, T)'
  matrix Xa = C*Xq
  matrix y = y_annual
  matrix u

  # initial OLS
  matrix b = mols(y, Xa, &u)
  scalar ac1 = corrgm(u,1)[1]
  matrix ma = {.5}
  scalar crit = BFGSmax(ma, chow_lin(&ma, ac1))
  matrix V = autocorr_matrix(ma[1], T4)
  matrix W = inv(qform(C,V))
  Xa = Xa'

  # GLS
  matrix beta = inv(qform(Xa,W))*Xa*W*y
  matrix uf = y - Xa'beta

  # resulting quarterly series
  matrix yq = Xq*beta + V*C'*W*uf
  return yq
end function

# expand annual series to quarterly by replicating
# the values for each quarter

function matrix expand_series (series y_annual, scalar T)
  matrix ystep = zeros(T*4, 1)
  scalar k = 1
  loop i=1..T -q
    loop j=1..4 -q
      ystep[k] = y_annual[i]
      k++
    endloop
  endloop
  return ystep
end function

# main script starts here

# open your data file
open data3-6
# and select your series for processing
series y = Ct

scalar T = $nobs
scalar T4 = T * 4
matrix yq = interpol_aq(y, T)

# for comparison, form simple expanded version of the
# original annual series
matrix ystep = expand_series(y, T)

# create quarterly dataset
nulldata T4 --preserve
setobs 4 1959:1
series C1 = 4*yq
series C2 = ystep
gnuplot C1 C2 --time-series --with-lines --output=display