set verbose off

function scalar beta1_SSR (scalar th2, const series y,
                           const series x, list L)
  matrix theta = {1, th2}
  series mdx = mlincomb(L, theta, 2)
  # run OLS conditional on theta
  ols y 0 x mdx --quiet
  return $ess
end function

function matrix midas_GNR (const matrix theta, const series y,
                           const series x, list L, int type)
  # Gauss-Newton regression
  series mdx = mlincomb(L, theta, type)
  ols y 0 x mdx --quiet
  matrix b = $coeff
  matrix u = {$uhat}
  matrix mgrad = mgradient(nelem(L), theta, type)
  matrix M = {const, x, mdx} ~ (b[3] * {L} * mgrad)
  matrix V
  set svd on # in case of strong collinearity
  mols(u, M, null, &V)
  return (b | theta) ~ sqrt(diag(V))
end function

/* main */

open gdp_midas.gdt --quiet

series dy = 100 * ldiff(qgdp)
series dy1 = dy(-1)
list dX = ld_payem*
list dXL = hflags(3, 11, dX)

# estimation sample
smpl 1985:1 2009:1

matrix b = {0, 1.01, 100}
# use Golden Section minimizer
SSR = GSSmin(b, beta1_SSR(b[1], dy, dy1, dXL), 1.0e-6)
printf "SSR (GSS) = %.15g\n", SSR
matrix theta = {1, b[1]}' # column vector needed
matrix bse = midas_GNR(theta, dy, dy1, dXL, 2)
bse[4,2] = $nan # mask std error of clamped coefficient
modprint bse "const dy(-1) HF_slope Beta1 Beta2"