How to apply seasonal decomposition and singular spectrum analysis (SSA) technique to time series

I am having a time series and want to apply seasonal decomposition technique to extract trend and seasonal component. After getting the seasonal component, I want to analyse the seasonal component using SSA.

Answer 1

[Pavl M.]

home
clear all
close all
pack


dt = readtable('time_series.txt');

%Prefilter - trim, + other datum munging if you'd like, which?
dt(isnan(dt)) = []; % Trims the NaN values from the vector
N = numel(dt) % Show the length of the new vector
lag = 5; %Which is suited for the model and decomposition requirements?
ind = dt.index;
d = dt.data;
[LT,ST,R] = trenddecomp(d,NumSeasonal=2); %SSA by default, you can select STL as well
figure
plot(ind,LT,'r',ST(:,1),'b',ST(:,2),'y',R,'g',d,'c')
legend('Long term','Seasonal','Residual','Original')

%Analysis of seasonal component using SSA?
%Correction found: SSA just extracts the seasonal components like STL, for
%anlysis you may want to use another strategies like cross-correlation,
%covariance... what else from statistics/machine learning to produce/serve
%for your which main objectives/goals?

%You may need additional SSA practice of:

% - creation of the trajectory matrix

M = 30;    % window length = embedding dimension

% - calculation of the covariance matrix

covX = xcorr(X,M-1,'unbiased');
Ctoep=toeplitz(covX(M:end));

% - eigendecomposition of the covariance matrix

Y=zeros(N-M+1,M);
for m=1:M
  Y(:,m) = X((1:N-M+1)+m-1);
end;
Cemb=Y'*Y / (N-M+1);

% - resulting eigenvalues, eigenvectors

C=Ctoep;

[RHO,LAMBDA] = eig(C);
LAMBDA = diag(LAMBDA);               % extract the diagonal elements
[LAMBDA,ind]=sort(LAMBDA,'descend'); % sort eigenvalues
RHO = RHO(:,ind);                    % and eigenvectors

figure
set(gcf,'name','Eigenvectors RHO and eigenvalues LAMBDA')
clf;
subplot(3,1,1);
plot(LAMBDA,'o-');
subplot(3,1,2);
plot(RHO(:,1:2), '-');
legend('1', '2');
subplot(3,1,3);
plot(RHO(:,3:4), '-');
legend('3', '4');

% - calculation of the principal components

PC = Y*RHO;

figure;
set(gcf,'name','Principal components PCs')
clf;
for m=1:4
  subplot(4,1,m);
  plot(t(1:N-M+1),PC(:,m),'k-');
  ylabel(sprintf('PC %d',m));
  ylim([-10 10]);
end;

% - reconstruction of the time series. 

RC=zeros(N,M);
for m=1:M
  buf=PC(:,m)*RHO(:,m)'; % invert projection
  buf=buf(end:-1:1,:);
  for n=1:N % anti-diagonal averaging
    RC(n,m)=mean( diag(buf,-(N-M+1)+n) );
  end
end;

figure
set(gcf,'name','Reconstructed components RCs')
clf;
for m=1:4
  subplot(4,1,m);
  plot(t,RC(:,m),'r-');
  ylabel(sprintf('RC %d',m));
  ylim([-1 1]);
end;

RC=zeros(N,M);
for m=1:M
  buf=PC(:,m)*RHO(:,m)'; % invert projection
  buf=buf(end:-1:1,:);
  for n=1:N % anti-diagonal averaging
    RC(n,m)=mean( diag(buf,-(N-M+1)+n) );
  end
end;

figure
set(gcf,'name','Reconstructed components RCs')
clf;
for m=1:4
  subplot(4,1,m);
  plot(t,RC(:,m),'r-');
  ylabel(sprintf('RC %d',m));
  ylim([-1 1]);
end;

%I did once cross-covariance and auto-covariance and Pearson and more
%correlative metrics, please ask me in private should customer will want to
%purchase it?