Constructing Multivariate Hodrick-Prescott Filters

Setup

Needs["EconGroup`EGPlotting`"];

Introduction

The standard (unvariate) HP filter finds a smoothed series based only on the time series properties of the original data.  It does so by finding the values of τ that minimise the objective function:

	ℒ = \(TraditionalForm
      \`∑\+\(t = 1\)\%T\((y\_t - τ\_t)\)\^2 + \ 
        λ 
          \(∑\+\(t = 2\)\%\(T - 1\)\((
                \((τ\_\(t + 1\) - τ\_t)\) - 
                  \((τ\_t - τ\_\(t - 1\))\))\)\^2\)\)

where the weight on smoothness (λ) is typically defined to be 1600 for quarterly data.  It turns out that the first-order conditions of this minimisation problem are linear equations that can be re-used for any series. They are also quite simple: other than the first two and last two conditions, they have the same form.  It is straightforward to write a function that finds HP-filtered trend of a series in Mathematica by using the general forms of these first-order conditions – the loss function does not have to be evaluated or differentiated each time, and the same function can be used for any value of λ.

\(HPFilter[data_?VectorQ, lambda_Integer] := 
	
      Module[{τ, focs, taulist, answer}, 
		
        With[{𝕋 = Length[data], λ = N[lambda]}, 
          
			taulist\  = \ 
            Table[τ[t], {t, 1, 𝕋}]; 

          focs = 	Thread[
              Flatten[Join[{
                      \(-2. \)\ 
                          data[[1]]\  + 
                        2. \ τ[1] + 2. \ λ\ τ[1] - 
                        4. \ λ\ τ[2] + 
                        2. \ λ\ τ[3], 
                      \(-2. \)\ 
                          data[[2]] - 
                        4. \ λ\ τ[1] + 2. \ τ[2] + 
                        10. \ λ\ τ[2] - 
                        8. \ λ\ τ[3] + 
                        2. \ λ\ τ[4]}, 

                    Table[	\(-2. \)\ \ 
                          data[[t]] + 
                        2. \ λ\ τ[t - 2] - 
                        8. \ λ\ τ[t - 1] + 2. \ τ[t] + 
                        12. \ λ\ τ[t] - 
                        8. \ λ\ τ[t + 1] + 
                        2. \ λ\ τ[t + 2], {t, 3, 
                        𝕋 - 2}], 

                    		{\(-2\)\ 
                          data[[𝕋 - 1
                              ]] + 
                        2. \ λ\ τ[𝕋 - 3] - 
                        8. \ λ\ τ[𝕋 - 2] + 
                        2. \ τ[𝕋 - 1] + 
                        10. \ λ\ 
                          τ[𝕋 - 1] - 
                        4. \ λ\ τ[𝕋], 
                      2. \ λ\ τ[𝕋 - 2] - 
                        4. \ λ\ τ[𝕋 - 1] + 
                        2. \ λ\ τ[𝕋]}]] == 
                0]; 
		answer = Solve[focs, taulist]; 
		
          Flatten[Expand[N[taulist /. answer]]]
					]]\)

Sometimes, however, a univariate approach isn’t enough: we may also want to use conditioning information from the series’ relationships with other variables to determine the true trend series. In that case, we need a multivariate Hodrick-Prescott filter. These are also quite simple to calculate using Mathematica: however, the structure of the first-order conditions will depend on the structure and number  of the conditioning equations, so there is no simple way to automate their determination.  Fortunately, the structure of the first-order conditions does not depend on the numerical parameter values in the conditioning equations, so you don’t have to re-write your HP filter functions each time you re-estimate the underlying equations and change the parameter values.

Calculating a Multivariate Hodrick-Prescott Filter

The loss function underlying a multivariate HP filter is a simple extension of the univariate case:

	ℒ = \(TraditionalForm
      \`∑\+\(i = 1\)\%n\( β\_i\) 
            \(∑\+\(t = 1\)\%T ϵ\_\(i, t\)\^2\) + 
        \(β\_trend\) 
          \(∑\+\(t = 1\)\%T\((y\_t - τ\_t)\)\^2\) + \ 
        λ 
          \(∑\+\(t = 2\)\%\(T - 1\)\((
                \((τ\_\(t + 1\) - τ\_t)\) - 
                  \((τ\_t - τ\_\(t - 1\))\))\)\^2\)\)

Where n is the number of conditioning equations,  the β’s are the weights on the deviations from trend and the conditioning equations, and the (Subscript[ϵ, i, t]) ’s are the residuals from the (i^th)conditioning equation at time t.

Finding the first-order conditions is as simple as defining a loss function like this, with algebraic symbols for the weights and the parameters in the conditioning equations, and differentiating it with respect to each of the (Subscript[τ, t]).  Because the first-order conditions of a Hodrick-Prescott filter – other than the first two and last two – have a common structure, we can use  a small horizon (T) in the loss function to find the correct structure for data series of any length.  The following section goes through this process step-by-step for a multivariate HP filter for output with three conditioning equations (prices, unemployment and capacity utilisation).

1. The Residuals

Firstly, we have to define the residuals, obtained by subtracting the body of the relavent behavioural equation from its dependent variable. This can be done with a standard Mathematica function definition. Note the [t_] construct instead of a subscript.  The “t_” translates as “anything, but for the purposes of the definition we’ll call it t”.  Also, this is all done in algebra: at this stage you don’t need data or numerical coefficients.  We also define the residual to be zero for values of  less than one: we assume that t=1 is the start of any series.

Definition for ϵπ

(Subscript[ϵ, p])[t_] := 0 /; t < 1(Subscript[ϵ, p])[t_] := 
      p[t] - (Subscript[δ, 1]) (y[t - 1] - τ[t - 1]) - 
        (Subscript[δ, 2]) (y[t - 2] - τ[t - 2])

Definition for ϵυ

(Subscript[ϵ, υ])[t_] := 0 /; t < 1(Subscript[ϵ, υ])[t_] :=  
      u[t] - θ (y[t - 1] - τ[t - 1]) - ξ

Definition for (Subscript[ϵ, c])

(Subscript[ϵ, c])[t_] := 0 /; t < 1(Subscript[ϵ, c])[t_] := 
      c[t]  -  (Subscript[α, 1]) (y[t - 1] - τ[t - 1]) - 
        (Subscript[α, 2]) (y[t - 2] - τ[t - 2])

Test these definitions

As these cells show, function definitions such as the ones above are “rewrite rules”: they tell Mathematica to substitute the right-hand side of the function definition whereever it encounters something matching the left-hand side of the function definition.

(Subscript[ϵ, p])[-3]

(Subscript[ϵ, c])[3]

2. The  Loss Function

The horizon over which to calculate the first-order conditions: 10 is usually enough, depending on the lag structure

𝕋 = 10;

Define the list of taus into a little vector

taulist = Table[τ[t], {t, 1, 𝕋}]

This is the definition of the loss function

You can see how Mathematica substitutes the definitions for the ϵ’s wherever they occur. Note that the terms in brackets are not expanded out automatically: you need to wrap an expression in the function ExpandAll[ ] for that to happen.

ℒ = 
      μ (∑_(t = 1)^𝕋) (Subscript[ϵ, p])[t]^2
               + β 
          (∑_(t = 1)^𝕋) (Subscript[ϵ, υ])[
                t]^2 + 
        ψ (∑_(t = 1)^𝕋) (Subscript[ϵ, c])[t]
                ^2 + 
        γ 
          (∑_(t = 1)^𝕋)(((
                y[t - 1] - τ[t - 1]))^2) + 
        λ 
          (∑_(t = 2)^(𝕋 - 1))(
              τ[t + 1]^2 + 4 τ[t]^2 + τ[t - 1]^2 - 
                4 τ[t + 1] τ[t] + 
                2 τ[t + 1] τ[t - 1] - 4 τ[t] τ[t - 1])

Finding the first-order conditions is as simple as telling Mathematica to take the partial derivative of the loss function above (ℒ) with respect to each of the τ[t]’s.  Note the use of the Expand[ ] function to expand out all the brackets in these derivatives. The first-order conditions are then just each of these derivatives equated to zero.

focs = 
      Table[{t, Expand[(∂_(τ[t]))ℒ]}, {t, 1, 
          𝕋}];

TableForm[focs]

You can see that these first-order conditions are still algebraic expressions.  That’s okay: because they are in algebraic form, you can use them as general definitions for any dataset.  With a little editing, they form the basis of the function we will construct to perform the filtering.

3. Constructing the function

This is the final function.  The structure is almost identical to the univariate HP filter function – not surprising since they are doing the same thing essentially, solving some linear equations given some data and parameter values.  The only differences are the structure of the first-order condition equations and the number of parameters that the user must supply (shown in the square brackets next to the function name at the beginning of the definition).

There are, however, a couple of catches to avoid when copying over the definitions for the first-order conditions into the Join statement in the function.

Firstly, function definitions don’t work if you have a variable with a subscript representing one of the parameters provided by the user: therefore, wherever (Subscript[α, 1]),  (Subscript[α, 2])   and the like appear in the derivatives above, something simpler (eg α1, α2) must be substituted.  Alternatively, one could use α1, α2 etc as the labels for those parameters in the original definitions of the residuals.

Secondly, the derivatives obtained above use single square brackets to denote time subscripts.  However, if you are dealing with vectors of numbers, you need to refer to elements of that vector (or Parts in Mathematica terminology), by using the part specification double-brackets instead of single brackets. For example Part i of vector X  is referred to like this:
	X[[i]], or if you want it to look pretty (and be easier to read) you can use the special double-brackets provided in Mathematica’s special mathematical fonts:   X[[i]].  Therefore, some editing of brackets to double-brackets will be required.

Finally, don’t forget that the Table[] statement covers the middle group of first-order conditions, so you need to re-edit the relevant first-order condition to be in a general form with the part specification in terms of the Table’s iterator variable (in this case, t).

\(multiHPFilter[y_?VectorQ, p_?VectorQ, u_?VectorQ, 
        c_?VectorQ, {λ_, μ_, ψ_, β_, γ_, 
          δ1_, δ2_, α1_, α2_, θ_, ξ_}]
         := 
	Module[{τ, focs, taulist, answer}, 
		
        With[{𝕋 = Length[y]}, 
			
          taulist\  = \ Table[τ[t], {t, 1, 𝕋}]; 

          focs = 	Thread[
              Flatten[Join[{
                      \(-2.\)\ β\ θ\ ξ + 
                        2.\ ψ\ 
                          c[[2]]\ 
                          α1 + 
                        2.\ ψ\ 
                          c[[3]]\ 
                          α2 + 
                        2.\ μ\ 
                          p[[2]]\ 
                          δ1 + 
                        2.\ μ\ 
                          p[[3]]\ 
                          δ2 + 
                        2.\ β\ θ\ 
                          u[[2]] - 
                        2.\ γ\ 
                          y[[1]] - 
                        2.\ β\ θ\^2\ 
                          y[[1]] - 
                        2.\ ψ\ α1\^2\ 
                          y[[1]] - 
                        2.\ ψ\ α2\^2\ 
                          y[[1]] - 
                        2.\ μ\ δ1\^2\ 
                          y[[1]] - 
                        2.\ μ\ \ δ2\^2\ 
                          y[[1]] - 
                        2.\ ψ\ α1\ α2\ 
                          y[[2]] - 
                        2.\ μ\ δ1\ δ2\ 
                          y[[2]] + 
                        2.\ γ\ τ[1] + 
                        2.\ β\ θ\^2\ τ[1] + 
                        2.\ λ\ τ[1] + 
                        2.\ ψ\ α1\^2\ τ[1] + 
                        2.\ ψ\ α2\^2\ τ[1] + 
                        2.\ μ\ δ1\^2\ τ[1] + 
                        2.\ μ\ δ2\^2\ τ[1] - 
                        4.\ λ\ τ[2] + 
                        2.\ ψ\ α1\ α2\ τ[2] + 
                        2.\ μ\ δ1\ δ2\ τ[2] + 
                        2.\ λ\ τ[3], 
						
                      \(-2\)\ β\ θ\ ξ + 	
                        2.\ ψ\ 
                          c[[3]]\ 
                          α1 + 
                        2.\ ψ\ 
                          c[[4]]\ 
                          α2 + 
                        2.\ μ\ 
                          p[[3]]\ 
                          δ1 + 
                        2.\ μ\ 
                          p[[4]]\ 
                          δ2 + 
                        2.\ β\ θ\ 
                          u[[3]] - 
                        2.\ ψ\ α1\ α2\ 
                          y[[1]] - 
                        2.\ μ\ δ1\ δ2\ 
                          y[[1]] - 
                        2.\ γ\ 
                          y[[2]] - 
                        2.\ β\ θ\^2\ 
                          y[[2]] - 
                        2\ ψ\ α1\^2\ 
                          y[[2]] - 
                        2.\ ψ\ α2\^2\ 
                          y[[2]] - 
                        2.\ μ\ δ1\^2\ 
                          y[[2]] - 
                        2.\ μ\ δ2\^2\ 
                          y[[2]] - 
                        2.\ ψ\ α1\ α2\ 
                          y[[3]] - 
                        2.\ μ\ δ1\ δ2\ 
                          y[[3]] - 
                        4.\ λ\ τ[1] + 
                        2.\ ψ\ α1\ α2\ τ[1] + 
                        2.\ μ\ δ1\ δ2\ τ[1] + 
                        2.\ γ\ τ[2] + 
                        2.\ β\ θ\^2\ τ[2] + 
                        10.\ λ\ τ[2] + 
                        2.\ ψ\ α1\^2\ τ[2] + 
                        2.\ ψ\ α2\^2\ τ[2] + 
                        2.\ μ\ δ1\^2\ τ[2] + 
                        2.\ μ\ δ2\^2\ τ[2] - 
                        8.\ λ\ τ[3] + 
                        2.\ ψ\ α1\ α2\ τ[3] + 
                        2.\ μ\ δ1\ δ2\ τ[3] + 
                        2.\ λ\ τ[4]}, 

                    Table[\(-2\)\ β\ θ\ ξ + 
                        2.\ ψ\ 
                          c[[t + 1]]\ 
                          α1 + 
                        2.\ ψ\ 
                          c[[t + 2]]\ 
                          α2 + 
                        2.\ μ\ 
                          p[[t + 1]]\ 
                          δ1 + 
                        2.\ μ\ 
                          p[[t + 2]]\ 
                          δ2 + 
                        2.\ β\ θ\ 
                          u[[t + 1]] - 
                        2.\ ψ\ α1\ α2\ 
                          y[[t - 1]] - 
                        2.\ μ\ δ1\ δ2\ 
                          y[[t - 1]] - 
                        2.\ γ\ 
                          y[[t]] - 
                        2.\ β\ θ\^2\ 
                          y[[t]] - 
                        2.\ ψ\ α1\^2\ 
                          y[[t]] - 
                        2.\ ψ\ α2\^2\ 
                          y[[t]] - 
                        2.\ μ\ δ1\^2\ 
                          y[[t]] - 
                        2.\ μ\ δ2\^2\ 
                          y[[t]] - 
                        2.\ ψ\ α1\ α2\^2\ 
                          y[[t + 1]] - 
                        2.\ μ\ δ1\ δ2\ 
                          y[[t + 1]] + 
                        2.\ λ\ τ[t - 2] - 
                        8.\ λ\ τ[t - 1] + 
                        2.\ ψ\ α1\ α2\ τ[t - 1] + 
                        2.\ μ\ δ1\ δ2\ τ[t - 1] + 
                        2.\ γ\ τ[t] + 
                        2.\ β\ θ\^2\ τ[t] + 
                        12.\ λ\ τ[t] + 
                        2.\ ψ\ α1\^2\ τ[t] + 
                        2.\ ψ\ α2\^2\ τ[t] + 
                        2.\ μ\ δ1\^2\ τ[t] + 
                        2.\ μ\ δ2\^2\ τ[t] - 
                        8.\ λ\ τ[t + 1] + 
                        2.\ ψ\ α1\ α2\ τ[t + 1] + 
                        2.\ μ\ δ1\ δ2\ τ[t + 1] + 
                        2.\ λ\ τ[t + 2], {t, 3, 
                        𝕋 - 2}], 

                    		{\(-2\)\ β\ θ\ ξ + 
                        2.\ ψ\ 
                          c[[𝕋]]\ α1 + 
                        2.\ μ\ 
                          p[[𝕋]] δ1 + 
                        2.\ β\ θ\ 
                          u[[𝕋]] - 
                        2.\ ψ\ α1\ α2\ 
                          y[[𝕋 - 2
                              ]] - 
                        2.\ μ\ δ1\ δ2\ 
                          y[[𝕋 - 2
                              ]] - 
                        2.\ γ\ 
                          y[[𝕋 - 1
                              ]] - 
                        2.\ β\ θ\^2\ 
                          y[[𝕋 - 1
                              ]] - 
                        2.\ ψ\ α1\^2\ 
                          y[[𝕋 - 1
                              ]] - 
                        2.\ μ\ δ1\^2\ 
                          y[[𝕋 - 1
                              ]] + 
                        2.\ λ\ τ[𝕋 - 3] - 
                        8.\ λ\ τ[𝕋 - 2] + 
                        2.\ ψ\ α1\ α2\ 
                          τ[𝕋 - 2] + 
                        2.\ μ\ δ1\ δ2\ 
                          τ[𝕋 - 2] + 
                        2.\ γ\ τ[𝕋 - 1] + 
                        2.\ β\ θ\^2\ 
                          τ[𝕋 - 1] + 
                        10.\ λ\ τ[𝕋 - 1] + 
                        2.\ ψ\ α1\^2\ 
                          τ[𝕋 - 1] + 
                        2.\ μ\ δ1\^2\ 
                          τ[𝕋 - 1] - 
                        4.\ λ\ τ[𝕋], 
                      2.\ λ\ τ[𝕋 - 2] - 
                        4.\ λ\ τ[𝕋 - 1] + 
                        2.\ λ\ τ[𝕋]}]] == 
                0]; 
\ \ \ \ \ \ \ \ answer = Solve[focs, taulist]; 
		
          Flatten[Expand[N[taulist /. answer]]]
					]]\)

4. A numerical example

A very important caveat!

A multivariate HP filter requires the user to choose the loss-function weights on the squared residuals from the conditioning equations (ie (Subscript[ϵ, p]),  (Subscript[ϵ, u]),  and so on).  It turns out that this is a very important decision: the data, coefficients and weights needs to be scaled so that one series does not dominate all other sources of variance.

Reading in the data

SetDirectory["HD256:"]

thedata = 
      ReadList["NaveenHPdata.csv", Word, 
        WordSeparators -> {","}];

data = 
      ToExpression /@Drop[Transpose[Drop[Partition[thedata, 5], 1]], 1];

Using the function

filtered  = 
      multiHPFilter[data[[1]], 
        data[[2]]/100., 
        data[[3]]/100., 
        data[[4]]/100., {1600., 100., 
          10., 10., 1., -0.0081, 0.1140, 0.0225, -0.0144058, 
          -0.72, .01105}]

A graph of the original and filtered series

Block[{graph1, graph2}, 
      graph1 = econplot[data[[1]], 
          "Multivariate HP Filter of output", 
          DisplayFunction -> Identity]; 
      graph2 = econplot[filtered, "Multivariate HP Filter of output", 
          PlotStyle -> {GrayLevel[0.5], AbsoluteThickness[2]}, 
          DisplayFunction -> Identity]; 
      Show[graph1, graph2, DisplayFunction -> $DisplayFunction]]