-
Notifications
You must be signed in to change notification settings - Fork 0
/
MW_test.m
59 lines (50 loc) · 1.68 KB
/
MW_test.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
function [Result,Err1,Err2] = MW_test(u,windowLength,threshold1,threshold2)
%
% [Result,Err1,Err2] =
% StationaryTestO2(u,windowLength,threshold1,threshold2,varargin) assess
% the first and second-order statioanrity of a random process.
%% Syntax
%
% Input:
%
% - u: double array [1xN] or [Nx1] : random process (time series) to test
%
% - windowLength: double [1x1]: Length of the moving window in time-steps (dimensionless)
%
% - threshold1: double [1x1]: Maximal acceptable value for the relative error between the
% isntantaneous mean and the static mean. Any value above threshold1
% will classify the time series as non-stationary.
%
% - threshold2: double [1x1]: Maximal acceptable value for the relative error between the
% isntantaneous standard deviation (std) and the static std.
% Any value above threshold2 will classify the time series as non-stationary.
%
% Output
%
% - Result: boolean: 1 if data is stationary. It is 0 otherwise
%
% - Err1: Calculated maximal relative error between the instantaneous mean and static
% mean
%
% - Err2: Calculated maximal relative error between the instantaneous std and static
% std
%
% Author: E. Cheynet - UiB - Last modified: 12-06-2020
%
[N1,N2]=size(u);
if N1>1 && N2>1
error('u must be a vector')
end
%% Moving mean function
movMeanU = movmean(u,windowLength);
Err1 = max(abs(movMeanU./nanmean(u)-1));
%% Moving standard deviation function
movStdU = movstd(detrend(u),windowLength); % I chose to remove the liner trend here
Err2 = max(abs(movStdU./nanstd(u)-1));
%% Assess the relative errors
if any(abs(Err1)>threshold1 | abs(Err2)>threshold2)
Result=0;
else
Result=1;
end
end