-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathExample1_WaSP_forecast.m
149 lines (123 loc) · 4.25 KB
/
Example1_WaSP_forecast.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
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
clc; clear
close all
% Created by Ze Jiang on 15/09/2021 ([email protected])
% Y: response = N x 1
% X: predictor= (N+N_fc) x n_var
N = 128; % number of observation
N_fc=24; % number of forecast (optional)
n_var=4; % number of variable
% initialize random number generator
%iseed = 101;
%rng(iseed,'twister')
%% synthetic data generation
fs = 50;
dt = 1/fs;
t = 0:dt:dt*(N+N_fc-1);
Y_ALL = (sin(2*pi*t+randn(1,1)) + 0.1*randn(size(t)))'; %+ t'; % sine wave add noise and trend
X = randn(N+N_fc,n_var);% + repmat(linspace(-pi,pi,N+N_fc)',1, n_var); % random predictors and trend
% X = [randn(N,n_var); repmat(100,N_fc,n_var)]; % test on abnormal value from forecast
Y = Y_ALL(1:N); Y_val = Y_ALL(N+1:end);
%% Daubechies wavelet with N vanishing moments
% N is a positive integer from 1 to 45
n_vanish = 1;
wname = ['db' num2str(n_vanish)] % db1 is equivalent to haar; 'sym2'; 'bior2.6'
flag_sign = 1;
%% method of discrete wavelet transform
switch 2
case 1
method = 'dwtmra' %dwtmra cannot be used in the forecast setting
case 2
method = 'modwt'
case 3
method = 'modwtmra'
case 4
method = 'at'
end
% maximum decomposition level: floor(log2(size(X,1)))
% or rule of thumb decomposition level: ceiling(log(n/(2*v-1))/log(2))-1 (Kaiser, 1994)
lev = floor(log2(size(X,1)))-1
%% WASP transform X_cal and forecasts together
% variance transformation using WaSP
[X_WaSP, C] = WaSP(Y, X, method, wname, lev, flag_sign);
% linear regression for each variable
RMSE=nan(1,n_var);
RMSE_WaSP=nan(1,n_var);
RMSE_opti=nan(1,n_var);
for i_var = 1:n_var
% optimal RMSE - Eq. 12 in WRR2020 paper
ratio=var(X(1:N,i_var))/var(X_WaSP(1:N,i_var));
%disp(ratio)
RMSE_opti(i_var) = sqrt((N-1)/N*(var(Y)-(norm(C(:,i_var))^2)*ratio));
% Std model: linear regression using original predictor
p_coeff1 = polyfit(X(1:N,i_var), Y, 1) ;
Y_fit = polyval(p_coeff1, X(1:N, i_var)) ;
RMSE(i_var) = sqrt(mean((Y-Y_fit).^2));
% VT model: linear regression using transformed predictor
p_coeff2 = polyfit(X_WaSP(1:N,i_var), Y, 1) ;
Y_fit = polyval(p_coeff2, X_WaSP(1:N, i_var)) ;
RMSE_WaSP(i_var) = sqrt(mean((Y-Y_fit).^2));
end
RMSE_WaSP
RMSE_opti
%% plot RMSE
figure
bar([RMSE',RMSE_WaSP',RMSE_opti']);
xlabel('No. of variable'); ylabel('RMSE');
legend('Std','VT','Optimal','NumColumns',1,'location','eastoutside')
title(['Variance transformation based on ',num2str(method) ' using ' num2str(wname)])
saveas(gca,'RMSE.fig');
%% plot - Y, X and X_WaSP
figure
sgtitle(['Calibration: ',num2str(method) ' using ' num2str(wname)])
for i_var = 1:n_var
subplot(n_var,1,i_var)
plot(Y, 'k')
xlim([0,N+N_fc]); %ylim([min(Y)*1.1,max(Y)*1.1]);
hold on
plot(X(:,i_var), 'b')
hold on
plot(X_WaSP(:,i_var), 'r')
hold off
legend('Y','X','Xnew','NumColumns',1,'location','eastoutside')
end
saveas(gca,'comparision.fig');
%% validation in the context of forecast
X_cal = X(1:N,:);
X_val = X(N+1:end,:); % predictors for validation
X_WaSP_cal = X_WaSP(1:N,:);
X_WaSP_val = X_WaSP(N+1:end,:);
% variance transformation using derived C from calibration period
if N_fc~=0
figure
sgtitle(['Validation: ',num2str(method) ' using ' num2str(wname)])
for i_var = 1:n_var
subplot(n_var,1,i_var)
plot(Y_val, 'k')
xlim([0,N_fc+1]);
hold on
plot(X_val(:,i_var), 'b')
hold on
plot(X_WaSP_val(:,i_var), 'r')
hold off
legend('Y','X','Xnew','NumColumns',1,'location','eastoutside')
end
saveas(gca,'comparision_val.fig');
% prediction
R_raw_fct = nan(1,n_var); R_WaSP_fct = nan(1,n_var);
R_WaSP_fct1 = nan(1,n_var);
for i_var = 1:n_var
% Std model: linear regression using original predictor
Y_fit1 = polyval(p_coeff1, X_val(:,i_var)) ;
R_raw_fct(i_var) = corr(Y_val,Y_fit1)^2;
% VT model: linear regression using transformed predictor
Y_fit2 = polyval(p_coeff2, X_WaSP_val(:,i_var)) ;
R_WaSP_fct(i_var) = corr(Y_val,Y_fit2)^2;
end
R2_mat = [R_raw_fct', R_WaSP_fct'];
figure
bar(R2_mat);
xlabel('No. of variable'); ylabel('R2');
legend('Std','WASP','NumColumns',1,'location','eastoutside')
title(['Variance transformation based on ',num2str(method) ' using ' num2str(wname)])
saveas(gca,'R2_val.fig');
end