hpfilterFS

hpfilterFS HP filter with missing/excluded observations via selection matrix.

expand all in page

Syntax

Description

Model interpretation (Gaussian).

$y=(y_1, y_2, \ldots, y_T)'$ is the observed time series. $y_bsb$ is a subset of units from $y$ of length $n_{bsb}$.

$Sy=y_{bsb} = Sm + S\epsilon$, $\epsilon \sim N(0, \sigma_\epsilon^2 I_{nbsb})$;

$D m = u$, $u \sim N(0, \sigma^2_\epsilon (1/\lambda) I_T)$.

HP ratio: $\lambda=\sigma^2_\epsilon/\sigma^2_u$, the greater $\lambda$, the smoother is the trend.

$D$ is the Second-difference matrix of size (T-2 x T):

each row has [1 -2 1].

$S$ is the selection matrix of size nbsb x T. Note that $SS'=I_{nbsb}$.

Conditioning on observed subset y_bsb = S y, using W = S'S (diag 0/1).

Posterior mean:

mhat = $\hat m= E(m|y_{bsb})=argmin_m ||S(y-m)||^2 + \lambda ||D m||^2$ = $(W + \lambda D'D) \ (W y)$.

Posterior cov:

$Cov(m|y_{bsb}) = \sigma^2_\epsilon (W + \lambda D'D)^{-1}$

example

out =hpfilterFS(y) Call to hpfilterFS with all default arguments.

example

out =hpfilterFS(y, Name, Value) Call to hpfilterFS with optional argument bsb.

Examples

expand all

  • Call to hpfilterFS with all default arguments.

    • Mdl = arima('Constant',0,'D',1,'MA',{0.5},'Variance',100);
n = 60;
y = simulate(Mdl,n);      
% In this case just the conditional mena 
out = hpfilterFS(y);

  • Call to hpfilterFS with optional argument bsb.

    • rng(1000)
Mdl = arima('Constant',0,'D',1,'MA',{0.5},'Variance',100);
n = 150;
y = simulate(Mdl,n); 
bsb = (1:round(n*0.9))';      
% In this case just the conditional mena 
out = hpfilterFS(y,'bsb',bsb);

    Related Examples

    expand all

  • Call to hpfilterFS with optional arguments bsb and predint.

    • rng(1000)
Mdl = arima('Constant',0,'D',1,'MA',{0.5},'Variance',100);
n = 150;
y = simulate(Mdl,n); 
bsb = (1:round(n*0.9))';      
% Prediction interval on all the observations (included and excluded)
out = hpfilterFS(y,'bsb',bsb,'predint','all','plots',true);
    Click here for the graphical output of this example (link to Ro.S.A. website)

  • Check equality with output of function hpfilter.

    • load Data_SchwertStock
TTM = rmmissing(DataTimeTableMth);
% Aggregate the monthly data in the timetable to quarterly measurements.
TTQ = convert2quarterly(TTM);
% Apply the Hodrick-Prescott filter to all variables in the quarterly timetable. The default smoothing parameter value is 1600. Display the last few observed components.
TQTT = hpfilter(TTQ);
TQTTchk=TQTT;
for j=1:size(TQTT,2)
outj=hpfilterFS(TTQ{:,j},'plots',0);
TQTTchk{:,j}=outj.mhat;
end
maxdiff=max(abs(TQTT{:,:}-TQTTchk{:,:}),[],"all");
assert(maxdiff<1e-10,"output of hpfilterFS different from hpfilter")

    Input Arguments

    expand all

    y — Time series to analyze. Vector or timetable.

    A row or a column vector with T elements, which contains the time series. Note that y may contain missing values. If y is a timetable then times of the timetable are shown in the plot of y and fitted values.

    Data Types: double or timetable

    Name-Value Pair Arguments

    Specify optional comma-separated pairs of Name,Value arguments. Name is the argument name and Value is the corresponding value. Name must appear inside single quotes (' '). You can specify several name and value pair arguments in any order as Name1,Value1,...,NameN,ValueN.

    Example: 'bsb',1:round(T/0.7): T=length(y); , 'conflev',0.999 , 'lambda',1000; , 'predint',"all"; , 's',52; , 'sigma2_eps',20; , 'ig_a0',10; , 'ig_a0',10; , 'plots',true;

    bsb —Indices of observed points to condition on.vector of length>=3.

    Indices of observed points used in the conditioning.

    If bsb is empty the the first 80 per cent units are used

    Example: 'bsb',1:round(T/0.7): T=length(y);

    Data Types: double

    conflev —confidence level for the confidence bands.scalar.

    A number between 0 and 1 which defines the confidence level which is used to produce the bands. The default value of conflev is 0.99.

    Example: 'conflev',0.999

    Data Types: double

    lambda —HP smoothing parameter.numeric scalar | empty value.

    The greater is the value of lambda the greater is the degree of smoothness which is used to estimate the trend. The default value of lambda which is used depends on the periodicity of the series (see optional input parameter s).

    If lambda is empty (default) the value of lambda which is used depends on the Ravn–Uhlig scaling rule (adjust by the 4th power of the observation frequency ratio):

    \[ \lambda(s)=1600(s/4​)^4 \]

    where $s$ is the number of observations per year Therefore, lambda=1600 for quarterly data, lambda=129600 for monthly data...

    Example: 'lambda',1000;

    Data Types: double or empty.

    predint —prediciton interval.string | empty value.

    String which specifies for which units to compute the variance of the prediction interval, that is the variance of the estimated trend plus the variance of the estimated noise component. If predint="nobsb" variance of the prediction is computed just for the units not beloing to bsb. If predint is "all" variance of the prediciton interval is computed for all the units. If predint="" variance of prediction interval is not computed.

    Example: 'predint',"all";

    Data Types: string or empty.

    s —length of seasonal period.numeric scalar | empty value.

    For monthly data s=12 (default), for quarterly data s=4,

    Example: 's',52;

    Data Types: double or empty.

    sigma2_eps —method to use to estimate residual variance.string of scalar numeric value.

    If sigma2_eps is a string possible values are: "MLaug" "MAPjef" "MAPig" "dfREML".

    Given:

    $RSS=||y_{bsb} - S \hat m||^2$.

    and Qhat:

    $\hat Q =RSS + \lambda ||D \hat m||^2$;

    $K = n_{bsb} + (n-2)$;

    $MLaug= ML (augmented) = \hat Q /K$.

    MAPjef= MAP (maximum a posteriori estimate) with Jeffreys prior.

    $MAPjef=\hat Q / (K+2)$.

    MAPig = MAP (maximum a posteriori estimate) with inverse gamma prior with parameters a0 and b0.

    $MAPig= (b_0 + 0.5\hat Q)/(a_0 + 1 + 0.5K)$.

    dfREML = df-REML-like (marginal smoother likelihood).

    $dfREML=RSS / (nbsb - df(lambda))$.

    $df(lambda)=trace(H) = trace(S A^{-1} S') = trace(A^{-1} S' S) = trace(A^{-1} W)$.

    The estimate of the $trace(A^{-1} W)$ is via Hutchinson:

    $tr(M) ≈ (1/niter) \sum_{i=1}^{niter} z'_i M z_i$.

    $z_i$ are Rademacher random. $niter$ is fixed to 50.

    variables, equal to +1 or −1 with equal probability.

    On the other hand, if sigma2_eps is a numeric scalar it is possible to supply the prior value.

    Example: 'sigma2_eps',20;

    Data Types: string or scalar double.

    ig_a0 —prior value of a in inverse gamma prior.numeric scalar.

    Prior value of parameter a of the inverse gamma to be used in case sigma2_eps="MAPig".

    Example: 'ig_a0',10;

    Data Types: scalar double.

    ig_b0 —prior value of b in IG.numeric scalar.

    Prior value of parameter b of the inverse gamma to be used in case sigma2_eps="MAPig".

    Example: 'ig_a0',10;

    Data Types: scalar double.

    plots —plot on the screeen of HP trend estimate.boolean.

    If plots = true a plot with the real time series with fitted values and trend estimate will appear on the screen. This plot is tagged forecastTS.

    The confidence bands which are shown depend on the input option predint.

    The default value of plot is false, that is no plot is shown on the screen.

    Example: 'plots',true;

    Data Types: logical.

    Output Arguments

    expand all

    out — description Structure

    A structure containing the following fields

    Value Description
    mhat

    n x 1 posterior mean of y. Estimated HP trend
    bsb

    indexes of the units used in the fit.
    sigma2_ML

    estimate of residual variance using augmented likelihood approach.
    sigma2_Jef

    MAP estimate of residual variance using Jeffreys prior.
    sigma2_IG

    MAP estimate of residual variance based on inverse gamma prior.
    sigma2_dfREML

    MAP estimate of residual variance based on marginal smoother likelihood.
    predVar

    n x 1 predictive variance or missing. This is the estimated trend variance + observation noise variance. The observation noise variance depends on the value of input option sigma2_eps. out.predVar is a scalar missing if input option predint is empty.

    out.PredVar is populated just for the units not belonging to bsb (input option predint is "nobsb"), or for all the units (input option predint is "all")
    PI_low

    n x 1 predictive variance or missing. This is the lower band of the confidence interval. out.PI_low is a scalar missing if input option predint is empty.

    out.PI_low is populated just for the units not belonging to bsb (input option predint is "nobsb"), or for all the units (input option predint is "all")
    PI_high

    n x 1 predictive variance or missing. This is the upper band of the confidence interval. out.PI_high is a scalar missing if input option predint is empty.

    out.PI_high is populated just for the units not belonging to bsb (input option predint is "nobsb"), or for all the units (input option predint is "all")

    References

    Hutchinson, M. F. (1989), A stochastic estimator of the trace of the influence matrix for laplacian smoothing splines, "Communications in Statistics Simulation and Computation", Vol. 18, pp. 1059–1076.

    Hodrick, R. J., & Prescott, E. C. (1997), Postwar U.S. Business Cycles:

    An Empirical Investigation, "Journal of Money, Credit and Banking", Vol. 29, pp. 1–16.

    Ravn, M. O., & Uhlig, H. (2002), On Adjusting the Hodrick–Prescott Filter for the Frequency of Observations, "Review of Economics and Statistics", Vol. 84, pp. 371–376.

    See Also

    | |

    This page has been automatically generated by our routine publishFS