Basic SSA (Singular Spectrum Analysis) implementation for Scilab.

ssa.sci

function [y_out, m, sig] = ssa_normalize(y)
	// http://en.wikipedia.org/wiki/Student%27s_t-statistic
	m = mean(y)
	sig = sqrt(variance(y))	// stdev(y)
	y_out = (y - m) / sig
endfunction
 
function [y_out] = ssa(y, L, I)
	[LAMBDA, U, V] = ssa_decompose(y, L)
	y_out = ssa_reconstruct(LAMBDA, U, V, I)
endfunction
 
function [LAMBDA, U, V, X] = ssa_decompose(y, L)
	// Ensure y is a column vector
	y = matrix(y, -1, 1)
 
// Stage 1: Decomposition
 
// Step 1: Embedding
 
	N = length(y)
	if L > N / 2 then
		L = N - L
	end
	K = N - L + 1
 
	// Time-delayed embedding of y, the trajectory matrix
	X = zeros(L, K)
	for i = 1 : K
		X(:, i) = y(i : i + L - 1)
	end
 
// Step 2: Singular value decomposition
 
	// Direct approach
	//[U, LAMBDA] = svd(X)
	//LAMBDA = diag(LAMBDA)
 
	// Covariance matrix
	C = X * X' / K
 
	// For stationary series
	//C = toeplitz(corr(X, L))
 
	[U, LAMBDA] = svd(C)
	// The eigenvalues of C are the squared eigenvalues of X
	LAMBDA = sqrt(diag(LAMBDA))
 
	// Principal components
	V = X' * U
	for i = 1 : L
		V(:, i) = V(:, i) / LAMBDA(i)
	end
endfunction
 
function [y] = ssa_reconstruct(LAMBDA, U, V, I)
	// Ensure I is a row vector
	I = matrix(I, 1, -1)
 
	[K, L] = size(V)
	N = K + L - 1
 
// Stage 2: Reconstruction
 
// Step 3: Grouping
 
	for i = I
		LAMBDA_U(:, i) = LAMBDA(i) * U(:, i)
	end
	// Reconstructed components
	//X = LAMBDA_U(:, I) * V_transposed(I, :)
	X = LAMBDA_U(:, I) * V(:, I)'
 
// Step 4: Diagonal averaging
 
	y = zeros(N, 1)
	for i = 1 : K + L - 1
		v = adiag(X, i)
		y(i) = sum(v) / length(v)
	end
	y = real(y)
endfunction
 
function [v] = adiag(x, z)
	// 2D matrix antidiagonal
 
	[N, M] = size(x)
 
	// The total diagonals count
	Z = N + M - 1
	// The current diagonal length
	z_len = min(z, Z - z + 1, min(N, M))
 
	i = max(z, N * (z - N + 1))
	step = N - 1
	v = x(i : step : i + step * z_len - 1)
endfunction