function series hp_via_kalman(series y, scalar lambda[0], bool oneside[0])
    if lambda == 0
        lambda = 100 * $pd^2
    endif
    # State transition matrix
    matrix T = {2, -1; 1, 0}
    # Observation matrix
    matrix Z = {1, 0}
    # Covariance matrix in the state equation
    matrix Q = {1/sqrt(lambda), 0; 0, 0}
    matrix my = {y}

    string desc = ""
    if oneside
        matrix my = my | 0
        desc = "1-sided "
    endif

    ssm = ksetup(my, Z', T, Q)
    ssm.obsvar = sqrt(lambda)
    ssm.inistate = {2*y[1]-y[2] ; 3*y[1]-2*y[2]}
    ssm.diffuse = 1

    err = oneside ? kfilter(&ssm) : ksmooth(&ssm)
    if err
        series ret = NA
    else
        mu = oneside ? ssm.state[2:,2] : ssm.state[,1]
        series ret = y - mu
    endif

    string d = sprintf("%sHP-filtered %s (lambda = %g)", desc, argname(y), lambda)
    setinfo ret --description="@d"
    return ret
end function

# --- example ------------------------------------------------------------

clear
open fedstl.bin
data houst
y = log(houst)

# one-sided, built-in then hansl
n1c = hpfilt(y, 1600, 1)
series h1c = hp_via_kalman(y, 1600, 1)
ols n1c const h1c --simple-print

# two-sided, built-in then hansl
n2c = hpfilt(y, 1600)
series h2c = hp_via_kalman(y, 1600)
ols n2c const h2c --simple-print