Back to the basics: don’t lag behind with your DataFrames.jl skills

By: Blog by Bogumił Kamiński

Re-posted from: https://bkamins.github.io/julialang/2022/09/23/lagging.html

Introduction

Several of my recent posts were a bit more advanced, so today I thought to write
something for beginners. A common task that one needs to do in data analysis
pipelines is lagging a variable in time. Today I will show you several patterns
how this can be done using DataFrames.jl and ShiftedArrays.jl. In particular, I
was prompted to write this post by a 2.0 release of the ShiftedArrays.jl package
and it reflects the current state of this package.

In the post I use Julia 1.8.1, DataFrames.jl 1.3.6, and ShiftedArrays.jl 2.0.

Setting up the scene

Let us load the packages that we will use in this post and create some starting
data frame with which we will work:

julia> using DataFrames

julia> using Dates

julia> using ShiftedArrays: lag

julia> df = DataFrame(day=Date.(2020, 9, [1, 2, 4, 5, 1, 2, 3, 6]),
                      id=repeat(1:2, inner=4),
                      v=1:8)
8×3 DataFrame
 Row │ day         id     v
     │ Date        Int64  Int64
─────┼──────────────────────────
   1 │ 2020-09-01      1      1
   2 │ 2020-09-02      1      2
   3 │ 2020-09-04      1      3
   4 │ 2020-09-05      1      4
   5 │ 2020-09-01      2      5
   6 │ 2020-09-02      2      6
   7 │ 2020-09-03      2      7
   8 │ 2020-09-06      2      8

Here an important pattern is using ShiftedArrays: lag. Since 2.0 release
the ShiftedArrays.jl package does not export the lag function. Therefore
to be able to use it without having to write ShiftedArrays.lag you need
to explicitly list it when loading the ShiftedArrays.jl package.

Simple lagging

Assume we want to lag variable v in the whole data frame. for this we can
use the transform function like this:

julia> transform(df, :v => lag)
8×4 DataFrame
 Row │ day         id     v      v_lag
     │ Date        Int64  Int64  Int64?
─────┼───────────────────────────────────
   1 │ 2020-09-01      1      1  missing
   2 │ 2020-09-02      1      2        1
   3 │ 2020-09-04      1      3        2
   4 │ 2020-09-05      1      4        3
   5 │ 2020-09-01      2      5        4
   6 │ 2020-09-02      2      6        5
   7 │ 2020-09-03      2      7        6
   8 │ 2020-09-06      2      8        7

If we wanted to lag it by more than one period we could write:

julia> transform(df, :v => (x -> lag(x, 2)) => :v_lag2)
8×4 DataFrame
 Row │ day         id     v      v_lag2
     │ Date        Int64  Int64  Int64?
─────┼───────────────────────────────────
   1 │ 2020-09-01      1      1  missing
   2 │ 2020-09-02      1      2  missing
   3 │ 2020-09-04      1      3        1
   4 │ 2020-09-05      1      4        2
   5 │ 2020-09-01      2      5        3
   6 │ 2020-09-02      2      6        4
   7 │ 2020-09-03      2      7        5
   8 │ 2020-09-06      2      8        6

As an exercise, here is an example how to do lagging by different periods
using a comprehension:

julia> transform(df, [:v => (x -> lag(x, i)) => "v_lag$i" for i in 1:3])
8×6 DataFrame
 Row │ day         id     v      v_lag1   v_lag2   v_lag3
     │ Date        Int64  Int64  Int64?   Int64?   Int64?
─────┼─────────────────────────────────────────────────────
   1 │ 2020-09-01      1      1  missing  missing  missing
   2 │ 2020-09-02      1      2        1  missing  missing
   3 │ 2020-09-04      1      3        2        1  missing
   4 │ 2020-09-05      1      4        3        2        1
   5 │ 2020-09-01      2      5        4        3        2
   6 │ 2020-09-02      2      6        5        4        3
   7 │ 2020-09-03      2      7        6        5        4
   8 │ 2020-09-06      2      8        7        6        5

Lagging by group

Simple lagging in our example seems not very useful. The reason is that it looks
that we should lag the :v variable in groups defined by the :id variable.
In DataFrames.jl this is simple, the only thing you need to change in the
above codes is passing groupby(df, :id) to transform:

julia> transform(groupby(df, :id), :v => lag)
8×4 DataFrame
 Row │ day         id     v      v_lag
     │ Date        Int64  Int64  Int64?
─────┼───────────────────────────────────
   1 │ 2020-09-01      1      1  missing
   2 │ 2020-09-02      1      2        1
   3 │ 2020-09-04      1      3        2
   4 │ 2020-09-05      1      4        3
   5 │ 2020-09-01      2      5  missing
   6 │ 2020-09-02      2      6        5
   7 │ 2020-09-03      2      7        6
   8 │ 2020-09-06      2      8        7

(the remaining examples can be updated using the same pattern)

Calendar lagging by group

The above example was better, but assumed that we want to lag data by
observations (so we shift a vector disregarding dates). What if we wanted to lag
the data by one calendar day. Here are two ways how you can do it (they slightly
differ in the produced result, so the choice of the approach depends on what you
would want in practice). In both approaches I assume that the passed dates are
unique by group.

In the first approach we want to keep all observations after lagging of the
:v variable and record only days that match some observations (before or after
lagging).

julia> combine(groupby(df, :id)) do sdf
           sort!(outerjoin(select(sdf, Not(:id)),
                           select(sdf,
                                  :day => ByRow(x -> x + Day(1)) => :day,
                                  :v => :v_lag),
                           on=:day), :day)
       end
12×4 DataFrame
 Row │ id     day         v        v_lag
     │ Int64  Date        Int64?   Int64?
─────┼─────────────────────────────────────
   1 │     1  2020-09-01        1  missing
   2 │     1  2020-09-02        2        1
   3 │     1  2020-09-03  missing        2
   4 │     1  2020-09-04        3  missing
   5 │     1  2020-09-05        4        3
   6 │     1  2020-09-06  missing        4
   7 │     2  2020-09-01        5  missing
   8 │     2  2020-09-02        6        5
   9 │     2  2020-09-03        7        6
  10 │     2  2020-09-04  missing        7
  11 │     2  2020-09-06        8  missing
  12 │     2  2020-09-07  missing        8

In this approach in each group defined by the :id column we create a helper
data frame with :day column incremented by one and outer join it with the
original data frame (with :id column dropped), next we sort the result on
:day.

In the second approach we keep only observations in the originally observed
period for a given group, but fill all the potential dates in the considered
period.

julia> combine(groupby(df, :id)) do sdf
           ref_df = DataFrame(day=minimum(sdf.day):Day(1):maximum(sdf.day))
           transform!(leftjoin!(ref_df, select(sdf, Not(:id)), on=:day),
                      :v => lag)
       end
11×4 DataFrame
 Row │ id     day         v        v_lag
     │ Int64  Date        Int64?   Int64?
─────┼─────────────────────────────────────
   1 │     1  2020-09-01        1  missing
   2 │     1  2020-09-02        2        1
   3 │     1  2020-09-03  missing        2
   4 │     1  2020-09-04        3  missing
   5 │     1  2020-09-05        4        3
   6 │     2  2020-09-01        5  missing
   7 │     2  2020-09-02        6        5
   8 │     2  2020-09-03        7        6
   9 │     2  2020-09-04  missing        7
  10 │     2  2020-09-05  missing  missing
  11 │     2  2020-09-06        8  missing

In this approach, again by :id group, we first create a reference data frame
with all dates in the observed period, next we leftjoin! the source data frame
into it (dropping :id as it is not needed), and finally just use the lag
function. Note that in this solution, in comparison to the previous one, rows
for dates 2020-09-06 (group 1) and 2020-09-07 (group 2) are missing, but we
have an extra row 2020-09-05 (group 2) filled with missing.

Conclusions

Lagging of variables is quite often needed in practice, especially in panel
econometrics context. I hope you will find the patterns I have presented in this
post useful in your data analysis projects. Additionally, I have chosen the
examples so that several functions from DataFrames.jl are employed in
combination.

Approximating a function using derivatives

By: Tamás K. Papp

Re-posted from: https://tamaspapp.eu/pages/blog/2022/hermite-approximation-spectralkit/index.html

When simulating from an economic model, I had to approximate a function \(f(x; \theta): [0,1] \to [0,1]\) for a variety of \(\theta\)s. \(f\) itself has to be solved for numerically, but otherwise it is pretty friendly, being continuous and increasing, with \(f(0)=0\) and \(f(1)=1\).

After profiling, this turned out to be the most costly part, so I had to approximate it. Since I needed derivatives \(f'(x)\), I was wondering whether making the approximation match them (known as a Hermite interpolation) would increase accuracy.

The (pedagogical, unoptimized) code below sums up the gist of my numerical experiments, with f below standing in for my implicitly solved function. It also demonstrates the new features of SpectralKit.jl v0.10.

First, we set up the problem:

using SpectralKit, PGFPlotsX, DisplayAsf(x) = (exp(x) - 1) / (exp(1) - 1)
f′(x) = exp(x) / (exp(1) - 1)
const I = BoundedLinear(0, 1)   # interval we map from

Then define an interpolation using N Chebyshev nodes, matching the values.

function interpolation0(f, N)
    basis = Chebyshev(EndpointGrid(), N)
    ϕ = collocation_matrix(basis) \ map(f ∘ from_pm1(I), grid(basis))
    linear_combination(basis, ϕ) ∘ to_pm1(I)
end;

Same exercise, but with the derivatives too, so we need two bases, one with double the number of functions (so we need to make sure N is even), while we just use N/2 for the nodes.

function interpolation01(f, f′, N)
    @assert iseven(N)
    basis1 = Chebyshev(EndpointGrid(), N ÷ 2) # nodes from this one
    basis2 = Chebyshev(EndpointGrid(), N)     # evaluate on this basis
    x = from_pm1.(I, grid(basis1))            # map nodes from [-1,1]
    M = collocation_matrix(basis2, to_pm1.(I, derivatives.(x)))
    ϕ = vcat(map(y -> y[0], M), map(y -> y[1], M)) \ vcat(f.(x), f′.(x))
    linear_combination(basis2, ϕ) ∘ to_pm1(I)
end;

Importantly, note that mapping to [-1,1] for the collocation matrix has to be preceded by lifting to derivatives.

Then calculate the max abs difference, in digits (log10).

function log10_max_abs_diff(f, f̂; M = 1000)
    x = range(0, 1; length = M)
    log10(maximum(@. abs(f(x) - f̂(x))))
end;

Then let's explore the errors in values …

Ns = 4:2:20
errors = [(log10_max_abs_diff(f, interpolation0(f, N)),
           log10_max_abs_diff(f, interpolation01(f, f′, N)))
          for N in Ns]
9-element Vector{Tuple{Float64, Float64}}:
 (-3.1996028783051695, -2.594513489315976)
 (-5.882145733021446, -5.488666848393999)
 (-8.835160643191552, -8.232779398084544)
 (-11.994023867372805, -11.44897859894945)
 (-15.176438519807359, -14.664555158828485)
 (-15.35252977886304, -15.35252977886304)
 (-15.255619765854984, -15.35252977886304)
 (-15.35252977886304, -15.35252977886304)
 (-15.35252977886304, -15.35252977886304)

… and derivatives.

d_errors = [(log10_max_abs_diff(f′, (x -> x[1]) ∘ interpolation0(f, N) ∘ derivatives),
             log10_max_abs_diff(f′, (x -> x[1]) ∘ interpolation01(f, f′, N) ∘ derivatives))
            for N in Ns]
9-element Vector{Tuple{Float64, Float64}}:
 (-2.0758500387125216, -2.093336352131656)
 (-4.549339116162139, -4.611253272379436)
 (-7.363367596306161, -7.429305371299876)
 (-10.417554370684012, -10.485171320264207)
 (-13.381718167990524, -13.689771947181466)
 (-13.834015838985154, -14.374806173574193)
 (-14.03551167781493, -14.539616422220185)
 (-13.724140848812729, -14.750469787535078)
 (-13.714040521908403, -14.724140848812729)

Finally the plots:

@pgf Axis({ xlabel = "number of basis functions",
            ylabel = "log10 abs error in values",
            legend_cell_align= "left" },
          PlotInc(Table(Ns, first.(errors))),
          LegendEntry("fitting values"),
          PlotInc(Table(Ns, last.(errors))),
          LegendEntry("fitting values and derivatives")) |> DisplayAs.SVG

@pgf Axis({ xlabel = "number of basis functions",
            ylabel = "log10 abs error in values",
            legend_cell_align= "left" },
          PlotInc(Table(Ns, first.(d_errors))),
          LegendEntry("fitting values"),
          PlotInc(Table(Ns, last.(d_errors))),
          LegendEntry("fitting values and derivatives")) |> DisplayAs.SVG

The conclusion is that even without matching them explicitly, derivatives are well-approximated. Getting an extra digit of accuracy in derivatives above 12–14 nodes means sacrificing a digit of accuracy with a low number of nodes. 14 seems to be the break-even point here, but then we are at machine precision anyway.

As usual, simply approximating with Chebyshev polynomials is extremely accurate in itself for practical purposes, even when derivatives are needed. Of course, this depends on the function being “nice”.

This page was generated using Literate.jl.
Download the source, project, and manifest.