Help on rolling window Fourier Transform and dydx

[Raymond Liu]

Hello there, I am new to the STATA Forum and I am hoping you could help me with generating a rolling window Fourier Transform and dydx. I have written a simple do-loop (see below) but it takes a while to generate if the sample size is very large. I wonder if there is a way to make it more efficient. I tried using "rolling" in STATA but it doesn't seem working.

clear
set obs 25600
gen t=(_n-1)*(2*_pi)/25600
gen y=exp(-(cos(t))^2)*(sin(2*t)+2*cos(4*t)+0.4*sin(t)*sin(50 *t))
egen seq = seq()
tsset seq
local x = 20
while `x' <= _N {
qui dydx y t if seq <= `x', gen(y1)
qui fft y if seq <= `x' & seq > `x' - 15, gen(v u)
if (`x' == 20) {
qui gen Slope_y = y1 if seq == `x'
qui gen fft1 = v if seq == `x'
}
else {
qui replace Slope_y = y1 if seq == `x'
qui replace fft1 = v if seq == `x'
}
qui drop y1 u v
local x = `x' + 1
dis `x'
}

Much appreciated.
Raymond

[Clyde Schechter]

The slowness is arising from the -if seq <= `x'- and -if seq <= `x' & seq > `x' - 15- clauses. Those conditions have to be re-checked for every observation in the entire data set on every iteration of the loop. You can replace these with equivalent -in- conditions that will be much, much faster.

First, the seq variable is entirely unnecessary. The way you have created it, seq will always have the same value as _n. This makes it possible to replace -if seq <= `x'- with -in 1/`x'-, and -if seq <= `x' & seq > `x'-15- with -in `=`x'-14'/`x'-, respectively. Making these changes will greatly speed up this code.

A few other tweaks might make barely noticeable improvements in speed if _N is really in the high millions:

1. Instead of

local `x' = 20
while `x' <= _N {
    ....
    local x = `x' + 1
}

use

forvalues x = 20/`=_N' {
    ....
}

2. Instead of distinguishing the `x' == 20 case with -gen-eration of y1 and fft1, vs -replace-ing them in the other cases, initialize both y1 and fft1 to missing values before the loop. Then -replace-ing them will work for both `x' == 20 and the other cases, eliminating a bit of code interpretation and branching.

[Clyde Schechter]

Hitting a slow patch in my day, I decided to try this on a nuts and bolts level. My statement that variable -seq- is unnecessary was wrong. You do need it so that you can -tsset seq- in preparation for -fft-. But you don't need -egen, seq()- for that: you can just -gen `c(obs_t)' seq = _n- for that. And then all the changes from -if- to -in- I proposed will still work. Also, in #2, I forgot to say that -replace ... if `x' == seq- should be changed to -replace ... in `x'-. The end result is:

clear
set obs 26500
gen t=(_n-1)*(2*_pi)/26500
gen y=exp(-(cos(t))^2)*(sin(2*t)+2*cos(4*t)+0.4*sin(t)*sin(50 *t))
gen `c(obs_t)' seq = _n
tsset seq
gen Slope_y = .
gen fft1 = .
forvalues x = 20/`=_N' {
   qui dydx y t in 1/`x', gen(y1)
   qui fft y in `=`x'-14'/`x', gen(v u)
    qui replace Slope_y = y1 in `x'
    qui replace fft1 = v in `x'
    qui drop y1 u  v
}

Nevertheless, the speedups obtained are very disappointing, far less than I imagined. Looking at the source code for -dydx-, I see that it uses -if- clauses internally, thereby undoing the efficiency I hoped would be gained by replacing those -if-'s with -in-'s. I don't know how fft works internally because it is a built-in command--but, it, too, might do scans of the entire data set to pick out its estimation sample, even when the code describes it with -in- instead of -if-. Even if it doesn't do that, I believe that the best FFT algorithms run in O(_N log _N) time.

Basically, any loop over observations that contains -if- clauses will run in O(_N²) time, whereas using -in- clauses instead makes it O(_N) time, because -if- is itself inherently O(_N) whereas -in- is O(1). But when, as here, the loop also contains other commands that themselves run in O(_N) time (or longer), then the loop will still take O(_N²) time (or longer).

So, sorry to say, I don't think there is much potential for a serious speedup here. I'm afraid I don't really have any good suggestions here.

[Raymond Liu]

Hi Clyde, Thank you so much for your speedy reply and helpful comment. This is fantastic! I did not realiise the potential efficiency difference between "if" and "in". I was thinking using Mata - but didn't find much improvement. I will give it a go.

Again, much appreciated.

Raymond

[Joseph Coveney]

. . . it takes a while to generate if the sample size is very large. I wonder if there is a way to make it more efficient.

If you profile your code (see below), you'll see that the overwhelming majority of the time spent in each iteration of the loop is with dydx. So speeding things up really needs to focus there.

If you look at the user's manual entry for dydx, you'll see that it generates a cubic spline for all observations and computes analytical derivatives from that—except for the first and last observations, for which it uses a quadratic approximation with the adjacent observations. Given that you throw away all of the derivatives except for the last observation's in each pass through the loop, you obviate the need for using dydx altogether and can just compute the quadratic approximation directly as shown below.

In the code shown below (incorporating Clyde's excellent suggestions, by the way), I've first abbreviated the scope to one-tenth of your sample size (I'm impatient) in order to show the code profiling and to prove that the direct calculation produces the identical result as dydx. I've attached the complete do-file and log file (file-name suffix is "_pilot").

I've also attached corresponding do-file and log file for the full-scale run with the direct computation instead of unnecessary calls to dydx (they omit the "_pilot" file-name suffix). The full-scale run (25 600 observations) completes in less than ten seconds.

version 18.0

clear *

quietly set obs 2560 // 25600
generate double t = (_n - 1) * (2 * _pi) / _N

#delimit ;
generate double y =
    exp(-cos(t)^2) *
    (
        sin(2 * t) +
        2 * cos(4 * t) +
        0.4 * sin(t) * sin(50 * t)
    );
#delimit cr

generate `c(obs_t)' seq = _n
tsset seq

quietly generate double Slope_y = .
quietly generate double fft1 = .
quietly generate double dydx = .

timer clear
forvalues x = 20/`=_N' {

    timer on 1
    quietly dydx y t in 1/`x', double generate(y1)
    quietly replace Slope_y = y1 in `x'
    timer off 1

    timer on 2
    quietly fft y in `=`x'-14'/`x', gen(v u)
    quietly replace fft1 = v in `x'
    timer off 2

    timer on 3
    local x1 = `x' - 1
    local x2 = `x' - 2
    quietly replace dydx = (y[`x'] - y[`x1']) / (t[`x'] - t[`x1']) + ///
        (y[`x'] - y[`x2']) / (t[`x'] - t[`x2']) - ///
        (y[`x1'] - y[`x2']) / (t[`x1'] - t[`x2']) in `x'
    timer off 3
    
    drop y1 u  v
}

timer list

quietly generate double del = Slope_y - dydx
summarize del, meanonly
assert r(max) <= 1e-14

exit

[Joseph Coveney]

quietly generate double del = Slope_y - dydx

Make that

quietly generate double del = abs(Slope_y - dydx)

[Raymond Liu]

Thanks Joseph! I will give it a try to do a direct computation. That sounds promising.

Raymond