Time Series Decomposition and Forecasting

// This source code is subject to the terms of the Mozilla Public License 2.0 at https://mozilla.org/MPL/2.0/
// © capissimo

//@version=4
study("Time Series Decomposition and Forecasting", overlay=true)
// for details see https://www.datavedas.com/smoothing-techniques-and-time-series-decomposition/

// EXPERIMENTAL

// General functions
print(data, txt) => 
    var lbl = label.new(bar_index, data, txt, xloc.bar_index, yloc.price, color.new(color.white,100), label.style_label_left, color.new(color.black,0), size.normal, text.align_left)
    label.set_xy(lbl, bar_index, data), label.set_text(lbl, txt)

// Stat module functions
mean(v) => array.avg(v)

std(v)  => array.stdev(v)

corr(v, w) => 
    std_v = std(v), std_w = std(w)
    std_v > 0 and std_w > 0 ? array.covariance(v, w) / std_v / std_w : 0 

best_fit(v, w) => // alpha = slope, beta = intercept
    beta  = corr(v, w) * std(w) / std(v)
    alpha = mean(w) - beta * mean(v)
    [alpha, beta]

getInts(sze) =>                 // vector with integer ids
    vec = array.new_int(sze, 0)
    for i=0 to sze-1
        array.set(vec, i, i)
    vec

getFloats(sze, floats) =>       // vector with floats
    vec = array.new_float(sze)
    for i=sze-1 to 0
        array.set(vec, i, nz(floats[i]))
    vec

calcComponents(v, n) => // n = season
    sze = array.size(v)
    ma = array.new_float(sze, 0)     // moving average (ma)
    for i=0 to sze-n
        slc = array.slice(v, i, i+n)
        array.set(ma, i, array.avg(slc))
    array.unshift(ma, 0)
    array.unshift(ma, 0)
    array.pop(ma)
    array.pop(ma)
    
    cma = array.new_float(sze, 0)    // centered ma
    si  = array.new_float(sze, 0)    // seasonality/irregularity 
    idx  = array.new_float(n, 0)     // seasonality index and 
    des  = array.new_float(sze, 0)   // deseasoned sample
    for i=2 to sze-3
        array.set(cma, i, (array.get(ma, i) + array.get(ma, i+1)) / 2)
        array.set(si,  i, array.get(v, i) / array.get(cma, i))
        array.set(idx, i%n, array.get(idx, i%n) + array.get(si, i))
    for i=0 to n-1
        array.set(idx,  i, array.get(idx, i) / 3)
    for i=0 to sze-1
        array.set(des,  i, array.get(v, i) /  array.get(idx, i%n))
    
    
    [ma, cma, si, idx, des]

trend(ints, al, bt) =>
    sze = array.size(ints)
    trd = array.new_float(sze, 0)
    for i=0 to sze-1
        array.set(trd, i, array.get(ints, i) * bt + al)
    trd
    
forecast(trd, idx) =>
    sze = array.size(trd)
    n   = array.size(idx)
    fct = array.new_float(sze, 0)
    for i=0 to sze-1
        array.set(fct, i, array.get(trd, i) * array.get(idx, i%n))
    fct

get_sample() =>  // sample of data from the above paper
    samp = array.new_float(0)
    array.push(samp,50),  array.push(samp,61), array.push(samp,73),  array.push(samp,84),  array.push(samp,59)
    array.push(samp,72),  array.push(samp,89), array.push(samp,100), array.push(samp,72),  array.push(samp,93)
    array.push(samp,110), array.push(samp,126), array.push(samp,80), array.push(samp,100), array.push(samp,123), array.push(samp,140)
    samp

// Inputs
lkb = input(16, 'Number of data points', options=[8,12,16,20,24,28,32,36,40,44,48]) 
n   = input(4,  'Season length', minval=4, maxval=4)  
dbg = input(true)

// Calculations
//samp = get_sample()          // get ficticious data sample
samp = getFloats(lkb, close)   // get real data sample

[ma, cma, si, idx, des] = calcComponents(samp, n)  // calculate ma, centered ma, seasonality/irregularity, seasonality index and deseasoned sample

// get variable ‘t’ which will be nothing but a time code in the form of concecutive indexes
t = getInts(lkb)             // get indexes of data items

[al, bt] = best_fit(t, des)  // next take regression of sample, getting its alpha & beta 

t2 = getInts(lkb+n)          // extend lookback window of 16 by a period of 4

trd = trend(t2, al, bt)      // extract trend

// forecast for the upcoming 4-part season. We do this by using the seasonality and trend component and by multiplying them.
fct = forecast(trd, idx)     // finally forecast including for the next period of 4

sz1 = array.size(ma), sz2 = array.size(cma), sz3 = array.size(si)  // array sizes for debugging purposes
sz4 = array.size(idx)
sz5 = array.size(des), sz6 = array.size(trd), sz7 = array.size(fct)

// Rendering
// find which bar is hours away from the current time
milliseconds_in_nhours = 1000 * 60 * 60 * n  // millisecs * secs * min * hours

var line ln = line.new(x1=time, y1=0, x2=time, y2=0, xloc=xloc.bar_time, color=color.new(color.blue,0), width=4)
// for i=sz2-2 to 1
//     line.set_xloc(ln, time[i-1],  time[i], xloc.bar_time)
//     line.set_y1(ln, array.get(cma, i-1))
//     line.set_y2(ln, array.get(cma, i))

if barstate.islast and dbg
    print(ohlc4, 'MA/4           ' + tostring(sz1) + ':' +tostring(ma)  +'\n'+
                 'Centered MA/4  ' + tostring(sz2) + ':' +tostring(cma) +'\n'+
                 'Seas & Irreg   ' + tostring(sz3) + ':' +tostring(si)  +'\n'+
                 'Seas Idx       ' + tostring(sz4) + ':' +tostring(idx) +'\n'+ 
                 'Deseas         ' + tostring(sz5) + ':' +tostring(des) +'\n'+
                 'Regression a/b ' + tostring(al)  +'--' +tostring(bt)  +'\n'+
                 'Trend          ' + tostring(sz6) + ':' +tostring(trd) +'\n'+
                 'Forecast       ' + tostring(sz7) + ':' +tostring(fct))
             
if barstate.islast 
    ln := line.new(x1=time[lkb], y1=array.get(fct, lkb-1), x2=time+milliseconds_in_nhours, y2=array.get(fct, 0), xloc=xloc.bar_time, color=color.new(color.black,0), width=3)
    line.delete(ln[1])