9
$\begingroup$

Background

In a paper from Epstein (1991): On obtaining daily climatological values from monthly means, the formulation and an algorithm for calculating Fourier interpolation for periodical and even-spaced values are given.

In the paper, the goal is to obtain daily values from monthly means by interpolation.

In short, it is assumed that unknown daily values can be represented by the sum of harmonic components: $$ y(t) = a_{0} + \sum_{j}\left[a_{j}\,\cos(2\pi jt/12)+b_{j}\,\sin(2\pi jt/12)\right] $$ In the paper $t$ (time) is expressed in months.

After some derviation, it is shown that the terms can be calculated by: $$ \begin{align} a_{0} &= \sum_{T}Y_{T}/12 \\ a_{j} &= \left[ (\pi j/12)/\sin(\pi j/12)\right] \times \sum_{T}\left[Y_{T}\,\cos(2\pi jT/12)/6 \right]~~~~~~~j=1,\ldots, 5 \\ b_{j} &= \left[ (\pi j/12)/\sin(\pi j/12)\right] \times \sum_{T}\left[Y_{T}\,\sin(2\pi jT/12)/6 \right]~~~~~~~j=1,\ldots, 5 \\ a_{6} &= \left[ (\pi j/12)/\sin(\pi j/12)\right]\times \sum_{T}\left[Y_{T}\cos(\pi T)/12\right] \\ b_{6} &= 0 \end{align} $$ Where $Y_{T}$ denote the monthly means and $T$ the month.

Harzallah (1995) summarizes this aproach as follows: "The interpolation is carried out by adding zeros to the spectral coefficients of data and by performing an inverse Fourier transform to the resulting extended coefficients. The method is equivalent to applying a rectangular filter to Fourier coefficients."


Questions

My goal is to use the above methodology for interpolation of weekly means to obtain daily data (see my previous question). In summary, I have 835 weekly means of count data (see the example dataset at the bottom of the question). There are quite a few things that I don't understand before I can apply the approach outlined above:

  1. How would the formulas have to be changed for my situation (weekly instead of monthly values)?
  2. How could the time $t$ be expressed? I assumed $t/835$ (or $t/n$ with $n$ data points in general), is that correct?
  3. Why does the author calculate 7 terms (i.e. $0\leq j \leq 6$)? How many terms would I have to consider?
  4. I understand that the question can probably be solved by using a regression approach and using the predictions for interpolation (thanks to Nick). Still, some things are unclear to me: How many terms of harmonics should be included in the regression? And what period should I take? How can the regression be done to ensure that the weekly means are preserved (as I don't want an exact harmonic fit to the data)?

Using the regression approach (which is also explained in this paper), I managed to get an exact harmonic fit to the data (the $j$ in my example would run through $1, \ldots, 417$, so I fitted 417 terms). How can this approach be modified -$~$if possible$~$- to achieve the conservation of the weekly means? Maybe by applying correction factors to each regression term?

The plot of the exact harmonic fit is:

Exact harmonic fit

EDIT

Using the signal package and the interp1 function, here's what I've managed to do using the example data set from below (many thanks to @noumenal). I use q=7 as we have weekly data:

# Set up the time scale

daily.ts <- seq(from=as.Date("1995-01-01"), to=as.Date("2010-12-31"), by="day")

# Set up data frame 

ts.frame <- data.frame(daily.ts=daily.ts, wdayno=as.POSIXlt(daily.ts)$wday,
                       yearday = 1:5844,
                       no.influ.cases=NA)

# Add the data from the example dataset called "my.dat"

ts.frame$no.influ.cases[ts.frame$wdayno==3] <- my.dat$case

# Interpolation

case.interp1 <- interp1(x=ts.frame$yearday[!is.na(ts.frame$no.influ.case)],y=(ts.frame$no.influ.cases[!is.na(ts.frame$no.influ.case)]),xi=ts.frame$yearday, method = c("cubic"))

# Plot subset for better interpretation
par(bg="white", cex=1.2, las=1)
plot((ts.frame$no.influ.cases)~ts.frame$yearday, pch=20,
     col=grey(0.4),
     cex=1, las=1,xlim=c(0,400), xlab="Day", ylab="Influenza cases")
lines(case.interp1, col="steelblue", lwd=1)

Cubicinterpo

There are two issues here:

  1. The curve seem to fit "too good": it goes through every point
  2. The weekly means are not conserved

Example dataset

structure(list(date = structure(c(9134, 9141, 9148, 9155, 9162, 
9169, 9176, 9183, 9190, 9197, 9204, 9211, 9218, 9225, 9232, 9239, 
9246, 9253, 9260, 9267, 9274, 9281, 9288, 9295, 9302, 9309, 9316, 
9323, 9330, 9337, 9344, 9351, 9358, 9365, 9372, 9379, 9386, 9393, 
9400, 9407, 9414, 9421, 9428, 9435, 9442, 9449, 9456, 9463, 9470, 
9477, 9484, 9491, 9498, 9505, 9512, 9519, 9526, 9533, 9540, 9547, 
9554, 9561, 9568, 9575, 9582, 9589, 9596, 9603, 9610, 9617, 9624, 
9631, 9638, 9645, 9652, 9659, 9666, 9673, 9680, 9687, 9694, 9701, 
9708, 9715, 9722, 9729, 9736, 9743, 9750, 9757, 9764, 9771, 9778, 
9785, 9792, 9799, 9806, 9813, 9820, 9827, 9834, 9841, 9848, 9855, 
9862, 9869, 9876, 9883, 9890, 9897, 9904, 9911, 9918, 9925, 9932, 
9939, 9946, 9953, 9960, 9967, 9974, 9981, 9988, 9995, 10002, 
10009, 10016, 10023, 10030, 10037, 10044, 10051, 10058, 10065, 
10072, 10079, 10086, 10093, 10100, 10107, 10114, 10121, 10128, 
10135, 10142, 10149, 10156, 10163, 10170, 10177, 10184, 10191, 
10198, 10205, 10212, 10219, 10226, 10233, 10240, 10247, 10254, 
10261, 10268, 10275, 10282, 10289, 10296, 10303, 10310, 10317, 
10324, 10331, 10338, 10345, 10352, 10359, 10366, 10373, 10380, 
10387, 10394, 10401, 10408, 10415, 10422, 10429, 10436, 10443, 
10450, 10457, 10464, 10471, 10478, 10485, 10492, 10499, 10506, 
10513, 10520, 10527, 10534, 10541, 10548, 10555, 10562, 10569, 
10576, 10583, 10590, 10597, 10604, 10611, 10618, 10625, 10632, 
10639, 10646, 10653, 10660, 10667, 10674, 10681, 10688, 10695, 
10702, 10709, 10716, 10723, 10730, 10737, 10744, 10751, 10758, 
10765, 10772, 10779, 10786, 10793, 10800, 10807, 10814, 10821, 
10828, 10835, 10842, 10849, 10856, 10863, 10870, 10877, 10884, 
10891, 10898, 10905, 10912, 10919, 10926, 10933, 10940, 10947, 
10954, 10961, 10968, 10975, 10982, 10989, 10996, 11003, 11010, 
11017, 11024, 11031, 11038, 11045, 11052, 11059, 11066, 11073, 
11080, 11087, 11094, 11101, 11108, 11115, 11122, 11129, 11136, 
11143, 11150, 11157, 11164, 11171, 11178, 11185, 11192, 11199, 
11206, 11213, 11220, 11227, 11234, 11241, 11248, 11255, 11262, 
11269, 11276, 11283, 11290, 11297, 11304, 11311, 11318, 11325, 
11332, 11339, 11346, 11353, 11360, 11367, 11374, 11381, 11388, 
11395, 11402, 11409, 11416, 11423, 11430, 11437, 11444, 11451, 
11458, 11465, 11472, 11479, 11486, 11493, 11500, 11507, 11514, 
11521, 11528, 11535, 11542, 11549, 11556, 11563, 11570, 11577, 
11584, 11591, 11598, 11605, 11612, 11619, 11626, 11633, 11640, 
11647, 11654, 11661, 11668, 11675, 11682, 11689, 11696, 11703, 
11710, 11717, 11724, 11731, 11738, 11745, 11752, 11759, 11766, 
11773, 11780, 11787, 11794, 11801, 11808, 11815, 11822, 11829, 
11836, 11843, 11850, 11857, 11864, 11871, 11878, 11885, 11892, 
11899, 11906, 11913, 11920, 11927, 11934, 11941, 11948, 11955, 
11962, 11969, 11976, 11983, 11990, 11997, 12004, 12011, 12018, 
12025, 12032, 12039, 12046, 12053, 12060, 12067, 12074, 12081, 
12088, 12095, 12102, 12109, 12116, 12123, 12130, 12137, 12144, 
12151, 12158, 12165, 12172, 12179, 12186, 12193, 12200, 12207, 
12214, 12221, 12228, 12235, 12242, 12249, 12256, 12263, 12270, 
12277, 12284, 12291, 12298, 12305, 12312, 12319, 12326, 12333, 
12340, 12347, 12354, 12361, 12368, 12375, 12382, 12389, 12396, 
12403, 12410, 12417, 12424, 12431, 12438, 12445, 12452, 12459, 
12466, 12473, 12480, 12487, 12494, 12501, 12508, 12515, 12522, 
12529, 12536, 12543, 12550, 12557, 12564, 12571, 12578, 12585, 
12592, 12599, 12606, 12613, 12620, 12627, 12634, 12641, 12648, 
12655, 12662, 12669, 12676, 12683, 12690, 12697, 12704, 12711, 
12718, 12725, 12732, 12739, 12746, 12753, 12760, 12767, 12774, 
12781, 12788, 12795, 12802, 12809, 12816, 12823, 12830, 12837, 
12844, 12851, 12858, 12865, 12872, 12879, 12886, 12893, 12900, 
12907, 12914, 12921, 12928, 12935, 12942, 12949, 12956, 12963, 
12970, 12977, 12984, 12991, 12998, 13005, 13012, 13019, 13026, 
13033, 13040, 13047, 13054, 13061, 13068, 13075, 13082, 13089, 
13096, 13103, 13110, 13117, 13124, 13131, 13138, 13145, 13152, 
13159, 13166, 13173, 13180, 13187, 13194, 13201, 13208, 13215, 
13222, 13229, 13236, 13243, 13250, 13257, 13264, 13271, 13278, 
13285, 13292, 13299, 13306, 13313, 13320, 13327, 13334, 13341, 
13348, 13355, 13362, 13369, 13376, 13383, 13390, 13397, 13404, 
13411, 13418, 13425, 13432, 13439, 13446, 13453, 13460, 13467, 
13474, 13481, 13488, 13495, 13502, 13509, 13516, 13523, 13530, 
13537, 13544, 13551, 13558, 13565, 13572, 13579, 13586, 13593, 
13600, 13607, 13614, 13621, 13628, 13635, 13642, 13649, 13656, 
13663, 13670, 13677, 13684, 13691, 13698, 13705, 13712, 13719, 
13726, 13733, 13740, 13747, 13754, 13761, 13768, 13775, 13782, 
13789, 13796, 13803, 13810, 13817, 13824, 13831, 13838, 13845, 
13852, 13859, 13866, 13873, 13880, 13887, 13894, 13901, 13908, 
13915, 13922, 13929, 13936, 13943, 13950, 13957, 13964, 13971, 
13978, 13985, 13992, 13999, 14006, 14013, 14020, 14027, 14034, 
14041, 14048, 14055, 14062, 14069, 14076, 14083, 14090, 14097, 
14104, 14111, 14118, 14125, 14132, 14139, 14146, 14153, 14160, 
14167, 14174, 14181, 14188, 14195, 14202, 14209, 14216, 14223, 
14230, 14237, 14244, 14251, 14258, 14265, 14272, 14279, 14286, 
14293, 14300, 14307, 14314, 14321, 14328, 14335, 14342, 14349, 
14356, 14363, 14370, 14377, 14384, 14391, 14398, 14405, 14412, 
14419, 14426, 14433, 14440, 14447, 14454, 14461, 14468, 14475, 
14482, 14489, 14496, 14503, 14510, 14517, 14524, 14531, 14538, 
14545, 14552, 14559, 14566, 14573, 14580, 14587, 14594, 14601, 
14608, 14615, 14622, 14629, 14636, 14643, 14650, 14657, 14664, 
14671, 14678, 14685, 14692, 14699, 14706, 14713, 14720, 14727, 
14734, 14741, 14748, 14755, 14762, 14769, 14776, 14783, 14790, 
14797, 14804, 14811, 14818, 14825, 14832, 14839, 14846, 14853, 
14860, 14867, 14874, 14881, 14888, 14895, 14902, 14909, 14916, 
14923, 14930, 14937, 14944, 14951, 14958, 14965, 14972), class = "Date"), 
    cases = c(168L, 199L, 214L, 230L, 267L, 373L, 387L, 443L, 
    579L, 821L, 1229L, 1014L, 831L, 648L, 257L, 203L, 137L, 78L, 
    82L, 69L, 45L, 51L, 45L, 63L, 55L, 54L, 52L, 27L, 24L, 12L, 
    10L, 22L, 42L, 32L, 52L, 82L, 95L, 91L, 104L, 143L, 114L, 
    100L, 83L, 113L, 145L, 175L, 222L, 258L, 384L, 755L, 976L, 
    879L, 846L, 1004L, 801L, 799L, 680L, 530L, 410L, 302L, 288L, 
    234L, 269L, 245L, 240L, 176L, 188L, 128L, 96L, 59L, 63L, 
    44L, 52L, 39L, 50L, 36L, 40L, 48L, 32L, 39L, 28L, 29L, 16L, 
    20L, 25L, 25L, 48L, 57L, 76L, 117L, 107L, 91L, 90L, 83L, 
    76L, 86L, 104L, 101L, 116L, 120L, 185L, 290L, 537L, 485L, 
    561L, 1142L, 1213L, 1235L, 1085L, 1052L, 987L, 918L, 746L, 
    620L, 396L, 280L, 214L, 148L, 148L, 94L, 107L, 69L, 55L, 
    69L, 47L, 43L, 49L, 30L, 42L, 51L, 41L, 39L, 40L, 38L, 22L, 
    37L, 26L, 40L, 56L, 54L, 74L, 99L, 114L, 114L, 120L, 114L, 
    123L, 131L, 170L, 147L, 163L, 163L, 160L, 158L, 163L, 124L, 
    115L, 176L, 171L, 214L, 320L, 507L, 902L, 1190L, 1272L, 1282L, 
    1146L, 896L, 597L, 434L, 216L, 141L, 101L, 86L, 65L, 55L, 
    35L, 49L, 29L, 55L, 53L, 57L, 34L, 43L, 42L, 13L, 17L, 20L, 
    27L, 36L, 47L, 64L, 77L, 82L, 82L, 95L, 107L, 96L, 106L, 
    93L, 114L, 102L, 116L, 128L, 123L, 212L, 203L, 165L, 267L, 
    550L, 761L, 998L, 1308L, 1613L, 1704L, 1669L, 1296L, 975L, 
    600L, 337L, 259L, 145L, 91L, 70L, 79L, 63L, 58L, 51L, 53L, 
    39L, 49L, 33L, 47L, 56L, 32L, 43L, 47L, 19L, 32L, 18L, 34L, 
    39L, 63L, 57L, 55L, 69L, 76L, 103L, 99L, 108L, 131L, 113L, 
    106L, 122L, 138L, 136L, 175L, 207L, 324L, 499L, 985L, 1674L, 
    1753L, 1419L, 1105L, 821L, 466L, 274L, 180L, 143L, 82L, 101L, 
    72L, 55L, 71L, 50L, 33L, 26L, 25L, 27L, 21L, 24L, 24L, 20L, 
    18L, 18L, 25L, 23L, 13L, 10L, 16L, 9L, 12L, 16L, 25L, 31L, 
    36L, 40L, 36L, 47L, 32L, 46L, 75L, 63L, 49L, 90L, 83L, 101L, 
    78L, 79L, 98L, 131L, 83L, 122L, 179L, 334L, 544L, 656L, 718L, 
    570L, 323L, 220L, 194L, 125L, 95L, 77L, 46L, 42L, 29L, 35L, 
    21L, 29L, 16L, 14L, 19L, 15L, 19L, 18L, 21L, 10L, 14L, 7L, 
    7L, 5L, 9L, 14L, 11L, 18L, 22L, 39L, 36L, 46L, 44L, 37L, 
    30L, 39L, 37L, 45L, 71L, 59L, 57L, 80L, 68L, 88L, 72L, 74L, 
    208L, 357L, 621L, 839L, 964L, 835L, 735L, 651L, 400L, 292L, 
    198L, 85L, 64L, 41L, 40L, 23L, 18L, 14L, 22L, 9L, 19L, 8L, 
    14L, 12L, 15L, 14L, 4L, 6L, 7L, 7L, 8L, 13L, 10L, 19L, 17L, 
    20L, 22L, 40L, 37L, 45L, 34L, 26L, 35L, 67L, 49L, 77L, 82L, 
    80L, 104L, 88L, 49L, 73L, 113L, 142L, 152L, 206L, 293L, 513L, 
    657L, 919L, 930L, 793L, 603L, 323L, 202L, 112L, 55L, 31L, 
    27L, 15L, 15L, 6L, 13L, 21L, 10L, 11L, 9L, 8L, 11L, 7L, 5L, 
    1L, 4L, 7L, 2L, 6L, 12L, 14L, 21L, 29L, 32L, 26L, 22L, 44L, 
    39L, 47L, 44L, 93L, 145L, 289L, 456L, 685L, 548L, 687L, 773L, 
    575L, 355L, 248L, 179L, 129L, 122L, 103L, 72L, 72L, 36L, 
    26L, 31L, 12L, 14L, 14L, 14L, 7L, 8L, 2L, 7L, 8L, 9L, 26L, 
    10L, 13L, 13L, 5L, 5L, 3L, 6L, 1L, 10L, 6L, 7L, 17L, 12L, 
    21L, 32L, 29L, 18L, 22L, 24L, 38L, 52L, 53L, 73L, 49L, 52L, 
    70L, 77L, 95L, 135L, 163L, 303L, 473L, 823L, 1126L, 1052L, 
    794L, 459L, 314L, 252L, 111L, 55L, 35L, 14L, 30L, 21L, 16L, 
    9L, 11L, 6L, 6L, 8L, 9L, 9L, 10L, 15L, 15L, 11L, 6L, 3L, 
    8L, 4L, 7L, 7L, 13L, 10L, 23L, 24L, 36L, 25L, 34L, 37L, 46L, 
    39L, 37L, 55L, 65L, 54L, 60L, 82L, 55L, 53L, 61L, 52L, 75L, 
    92L, 121L, 170L, 199L, 231L, 259L, 331L, 357L, 262L, 154L, 
    77L, 34L, 41L, 21L, 17L, 16L, 7L, 15L, 11L, 7L, 5L, 6L, 13L, 
    7L, 6L, 8L, 7L, 1L, 11L, 9L, 3L, 9L, 9L, 8L, 15L, 19L, 16L, 
    10L, 12L, 26L, 35L, 35L, 41L, 34L, 30L, 36L, 43L, 23L, 55L, 
    107L, 141L, 217L, 381L, 736L, 782L, 663L, 398L, 182L, 137L, 
    79L, 28L, 26L, 16L, 14L, 8L, 4L, 4L, 6L, 6L, 11L, 4L, 5L, 
    7L, 7L, 6L, 8L, 2L, 3L, 3L, 1L, 1L, 3L, 3L, 2L, 8L, 8L, 11L, 
    10L, 11L, 8L, 24L, 25L, 25L, 33L, 36L, 51L, 61L, 74L, 92L, 
    89L, 123L, 402L, 602L, 524L, 494L, 406L, 344L, 329L, 225L, 
    136L, 136L, 84L, 55L, 55L, 42L, 19L, 28L, 8L, 7L, 2L, 7L, 
    6L, 4L, 3L, 5L, 3L, 3L, 0L, 1L, 2L, 3L, 2L, 1L, 2L, 2L, 9L, 
    4L, 9L, 10L, 18L, 15L, 13L, 12L, 10L, 19L, 15L, 22L, 23L, 
    34L, 43L, 53L, 47L, 57L, 328L, 552L, 787L, 736L, 578L, 374L, 
    228L, 161L, 121L, 96L, 58L, 50L, 37L, 14L, 9L, 6L, 15L, 12L, 
    9L, 1L, 6L, 4L, 7L, 7L, 3L, 6L, 9L, 15L, 22L, 28L, 34L, 62L, 
    54L, 75L, 65L, 58L, 57L, 60L, 37L, 47L, 60L, 89L, 90L, 193L, 
    364L, 553L, 543L, 676L, 550L, 403L, 252L, 140L, 125L, 99L, 
    63L, 63L, 76L, 85L, 68L, 67L, 38L, 25L, 24L, 11L, 9L, 9L, 
    4L, 8L, 4L, 6L, 5L, 2L, 6L, 4L, 4L, 1L, 5L, 4L, 1L, 2L, 2L, 
    2L, 2L, 3L, 4L, 4L, 7L, 5L, 2L, 10L, 11L, 17L, 11L, 16L, 
    15L, 11L, 12L, 21L, 20L, 25L, 46L, 51L, 90L, 123L)), .Names = c("date", 
"cases"), row.names = c(NA, -835L), class = "data.frame")
$\endgroup$
  • 2
    $\begingroup$ Older literature often delights in formulas of rococo complexity for this kind of thing. In practice I would lean over backwards to thinking of this as a regression problem so that the interpolated values are just the predicted values from the regression. See e.g. stats.stackexchange.com/questions/60500/… The key principle is that the main cycle repeats once per year. $\endgroup$ – Nick Cox Jul 3 '13 at 14:29
  • 1
    $\begingroup$ In practice you want a very close fit to the data because you want local smoothing. That might need several Fourier pairs, but diminishing returns set in very quickly, so that each new pair of terms soon adds very little. You have just have to suck it and see. Plotting everything makes it clearer. $\endgroup$ – Nick Cox Jul 3 '13 at 15:36
  • 1
    $\begingroup$ I tried this out briefly with your data (using Stata, not R). Briefly, although there is marked seasonality in your data it is not regular enough for this approach to work well at all: for example, not only does the timing of the peaks vary a lot, but also the number of cases at the peak; in some years but not all there is a secondary peak late in the calendar year. Further, the seasonality is compounded with marked long-term trend. My guess is that to get daily cases you should resort to some strictly local interpolation or smoothing of the weekly series expanded seven-fold. $\endgroup$ – Nick Cox Jul 5 '13 at 15:21
  • 1
    $\begingroup$ In control systems engineering the Nyquist criterion is used as a floor for the sampling rate. It says sample at greater than twice the highest frequency in your data. In practice it is more conventional to sample above 5x the highest frequency you hope to resolve. If your input is weekly data, then Nyquist suggests the highest resolvable frequency is on the order of 2 weeks. It might be better if you had other weekly statistics to inform the sampling and buttress the mean. en.wikipedia.org/wiki/Nyquist-Shannon_sampling_theorem $\endgroup$ – EngrStudent Jul 8 '13 at 12:24
  • $\begingroup$ +1 Excellent question! Do you know by any chance how to detect signal in the noise (preferably in R), provided you have a number of distributions, where noise is a Gaussian and signal is a Gaussian plus another part? I looked at many packages and functions (signal, fft(), etc.) and even played with data, trying to apply Fourier transformation and even entropy measures, but to no avail so far. I was trying to answer a question (not mine) and learn something new along the way, as I find the topic quite interesting. $\endgroup$ – Aleksandr Blekh Jan 31 '15 at 4:40
4
+25
$\begingroup$

I am no expert on Fourier transforms, but...

Epstein's total sample range was 24 months with a monthly sample rate: 1/12 years. Your sample range is 835 weeks. If your goal is to estimate the average for one year with data from ~16 years based on daily data you need a sample rate of 1/365 years. So substitute 52 for 12, but first standardize units and expand your 835 weeks to 835*7 = 5845 days. However, if you only have weekly data points I suggest a sample rate of 52 with a bit depth of 16 or 17 for peak analysis, alternatively 32 or 33 for even/odd comparison. So the default input options include: 1) to use the weekly means (or the median absolute deviation, MAD, or something to that extent) or 2) to use the daily values, which provide a higher resolution.

Liebman et al. chose the cut-off point jmax = 2. Hence, Fig 3. contains fewer partials and is thus more symmetrical at the top of the sine compared to Fig 2. (A single partial at the base frequency would result in a pure sine wave.) If Epstein would have selected a higher resolution (e.g. jmax = 12) the transform would presumably only yield minor fluctuations with the additional components, or perhaps he lacked the computational power.

Through visual inspection of your data you appear to have 16-17 peaks. I would suggest you set jmax or the "bit depth" to either 6, 11, 16 or 17 (see figure) and compare the outputs. The higher the peaks, the more they contribute to the original complex waveform. So assuming a 17-band resolution or bit depth the 17th partial contributes minimally to the original waveform pattern compared to the 6th peak. However, with a 34 band-resolution you would detect a difference between even and odd peaks as suggested by the fairly constant valleys. The bit depth depends on your research question, whether you are interested in the peaks only or in both peaks and valleys, but also how exactly you wish to approximate the original series.

The Fourier analysis reduces your data points. If you were to inverse the function at a certain bit depth using a Fourier transform you could probably cross-check if the new mean estimates correspond to your original means. So, to answer your fourth question: the regression parameters you mentioned depend on the sensitivity and resolution that you require. If you do not wish for an exact fit, then by all means simply input the weekly means in the transform. However, beware that lower bit depth also reduces the data. For example, note how Epstein's harmonic overlay on Lieberman and colleagues' analysis misses the mid-point of the step function, with a skewed curve slightly to the right (i.e. temp. est. too high), in December in Figure 3.


Liebman and Colleagues' Parameters:

  • Bit Depth: 2

Epstein's Parameters:

  • Sample Rate: 12 [every month]
  • Sample Range: 24 months
  • Bit Depth: 6

Your Parameters:

  • Sample Rate: 365 [every day]

  • Sample Range: 5845 days

Exact Bit Depth Approach

Exact fit based on visual inspection. (If you have the power, just see what happens compared to lower bit-depths.)

  • Full Spectrum (peaks): 17
  • Full Spectrum (even/odd): 34

Variable Bit Depth Approach

This is probably what you wish to do:

  • Compare Peaks Only: 6, 11, 16, 17
  • Compare Even/Odd: 12, 22, 32, 34
  • Resynthesize and compare means

This approach would yield something similar to the comparison of Figures in Epstein if you inverse the transformation again, i.e. synthesise the partials into an approximation of the original time series. You could also compare the discrete points of the resynthesized curves to the mean values, perhaps even test for significant differences to indicate the sensitivity of your bit depth choice.


UPDATE 1:

Bit Depth

A bit - short for binary digit - is either 0 or 1. The bits 010101 would describe a square wave. The bit depth is 1 bit. To describe a saw wave you would need more bits: 0123210. The more complex a wave gets the more bits you need:

This is a somewhat simplified explanation, but the more complex a time series is, the more bits are required to model it. Actually, "1" is a sine wave component and not a square wave (a square wave is more like 3 2 1 0 - see attached figure). 0 bits would be a flat line. Information gets lost with reduction of bit depth. For example, CD-quality audio is usually 16 bit, but land-line phone quality audio is often around 8 bits.

Please read this image from left to right, focusing on the graphs:

FFT

You have actually just completed a power spectrum analysis (although at high resolution in your figure). Your next goal would be to figure out: How many components do I need in the power spectrum in order to accurately capture the means of the time series?

UPDATE 2


To Filter or not to Filter

I am not entirely sure how you would introduce the constraint in the regression as I am only familiar with interval constraints, but perhaps DSP is your solution. This is what I figured so far:

  • Step 1. Break down the series into sinus components through Fourier function on the complete data set (in days)

  • Step 2. Recreate the time series through an inverse Fourier transform, with the additional mean-constraint coupled to the original data: the deviations of the interpolations from the original means should cancel out each other (Harzallah, 1995).

My best guess is that you would have to introduce autoregression if I understand Harzallah (1995, Fig 2) correctly. So that would probably correspond to an infinite response filter (IIR)?

IIR http://paulbourke.net/miscellaneous/ar/

In summary:

  1. Derive means from Raw data
  2. Fourier Transform Raw data
  3. Fourier Inverse Transform transformed Data.
  4. Filter the result using IIR

Perhaps you could use an IIR filter without going through the Fourier analysis? The only advantage of the Fourier analysis as I see it is to isolate and determine which patterns are influential and how often they do reoccur (i.e. oscillate). You could then decide to filter out the ones that contribute less, for example using a narrow notch filter at the least contributing peak (or filter based on your own criteria). For starters, you could filter out the less contributing odd valleys that appear more like noise in the "signal". Noise is characterized by very few cases and no pattern. A comb filter at odd frequency components could reduce the noise - unless you find a pattern there.

Here's some arbitrary binning—for explanatory purposes only: Can you see the noise in the valleys?

Oops - There's an R Function for that!?

When searching for an IIR-filter I happen to discovered the R functions interpolate in the Signal package. Forget everything I said up to this point. The interpolations should work like Harzallah's: http://cran.r-project.org/web/packages/signal/signal.pdf

Play around with the functions. Should do the trick.


UPDATE 3

interp1 not interp

case.interp1 <- interp1(x=(ts.frame$no.influ.cases[!is.na(ts.frame$no.influ.case)]),y=ts.frame$yearday[!is.na(ts.frame$no.influ.case)],xi=mean(WEEKLYMEANSTABLE),method = c("cubic"))

Set xi to the original weekly means.

$\endgroup$
  • $\begingroup$ Thank you very much for this answer! My research goal is simple: I have weekly means and want to get daily estimates and the mean of the (interpolated) daily estimates over one week should equal the weekly mean (i.e. the original data point). Do you think that's possible?. Further, I don't understand what "bit depth" and "peak analysis" mean (I have no experience with Fourier transformations whatsoever). $\endgroup$ – COOLSerdash Jul 11 '13 at 19:35
  • 1
    $\begingroup$ @COOLSerdash see my update. Yes it is possible! But you would need to figure out the best way to compare the estimated means from the resynthesized time series with the original means from the original time series. $\endgroup$ – noumenal Jul 11 '13 at 19:53
  • $\begingroup$ (Btw: +1 tomorrow, I can't vote today anymore). Thanks a lot for the update, it is clearer now. I thought about the following procedure: 1) fit a fourier function to the weekly means by regression, 2) use the predictions of the regression to "fill in" the gaps between the weekly values (i.e. to get the daily values) 3) for each week, calculate the mean of all daily values and this mean should equal the original value. In the paper, Epstein used some kind of "correction" factor to force the function to have the desired properties, but I'm still not sure how to do that with regression. $\endgroup$ – COOLSerdash Jul 11 '13 at 20:01
  • $\begingroup$ @COOLSerdash See update 2! Skip to the final paragraph. $\endgroup$ – noumenal Jul 11 '13 at 21:10
  • $\begingroup$ Absolutely fantastic! Thank you very much for the research. Note that I've managed to implement Harzallah's approach using splines (linear and cubic). So I guess I would need interp. I've edited my question. Thanks again very much. $\endgroup$ – COOLSerdash Jul 11 '13 at 22:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.