Seasonal decomposition of a daily time series

I have a time series that contains daily variation from 2015 to 2023. I want to see the trend, seasonal components and unmodelled part. Kindly help. The time series has veen attached

Answer 1

[Steven Lord]

Take a look at the trenddecomp function (introduced in release R2021b) and/or the Find and Remove Trends Live Editor Task (introduced in release R2019b.)

[Hfile]

these two commmands do not work

[Hfile]

my version is older

[Steven Lord]

Which release are you using?

You could try using the detrend function as part of writing your own seasonal decomposition code. That function was introduced sometime prior to release R2006a.

Answer 2

[Star Strider]

Your data are not regularly-sampled, so it is necessary to use the funciton on them first. After that, calculate the Fourier transform of the data to see the relevant peaks. From there, you can use the lowpass or bandpass functions as appropriate to get the relevant information from your data.

Try this —

format longG
A1 = readmatrix('time_series.txt.txt')

A1 = 2465x2

1.0e+00 * 2015.00067971 0.169734 2015.00341756 0.168658 2015.00615541 0.16273 2015.00889326 0.166784 2015.01163111 0.16894 2015.01436896 0.162698 2015.01710681 0.165437 2015.01984466 0.165807 2015.02258252 0.163822 2015.02532037 0.166965

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Ts = mean(diff(A1(:,1)))

Ts = 
       0.00364454163149355

Fs = 1/Ts

Fs = 
          274.382926884058

Tsd = std(diff(A1(:,1)))

Tsd = 
        0.0101315392133316

[A2(:,2), A2(:,1)] = resample(A1(:,2), A1(:,1), Fs);
A2

A2 = 2465x2

1.0e+00 * 2015.00067971 0.169734 2015.00432425163 0.166694828901561 2015.00796879326 0.165415119984054 2015.01161333489 0.168926002473624 2015.01525787653 0.163587289904377 2015.01890241816 0.165679663045187 2015.02254695979 0.163847781821886 2015.02619150142 0.1664769105871 2015.02983604305 0.165745285427315 2015.03348058468 0.16090356601093

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

figure
plot(A1(:,1), A1(:,2))
grid
xlabel('Time (years)')
ylabel('Amplitude')
title('Original')

figure
plot(A2(:,1), A2(:,2))
grid
xlabel('Time (years)')
ylabel('Amplitude')
title('Resampled')

[FTs1,Fv] = FFT1(A2(:,2), A2(:,1));

[pks,locsidx] = findpeaks(abs(FTs1)*2, 'MinPeakProminence',2E-3)

pks = 3x1

0.00472265164707644 0.0111256141963358 0.0022965835898268

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

locsidx = 3x1

3 16 30

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

figure
plot(Fv, abs(FTs1)*2)
hold on
plot(Fv(locsidx), pks, 'vr')
hold off
grid
xlim([0 10])
xlabel('Frequency (Years^{-1})')
ylabel('Magnitude (Units)')
text(Fv(locsidx), pks, compose('\\leftarrow Frequency = %.3f years^{-1} =  (1/%.3f years)', [Fv(locsidx); 1./Fv(locsidx)].') )

function [FTs1,Fv] = FFT1(s,t)
% Arguments:
%    s: Signal Vector Or Matrix
%    t: Associated Time Vector 
    t = t(:);
    L = numel(t);
    if size(s,2) == L
        s = s.';
    end
    Fs = 1/mean(diff(t));
    Fn = Fs/2;
    NFFT = 2^nextpow2(L);
    FTs = fft((s - mean(s)) .* hann(L).*ones(1,size(s,2)), NFFT)/sum(hann(L));
    Fv = Fs*(0:(NFFT/2))/NFFT;
%     Fv = linspace(0, 1, NFFT/2+1)*Fn;
    Iv = 1:numel(Fv);
    FTs1 = FTs(Iv,:);
end

Experiment to get the result you want. I will help with the filtering if necessary.

.