-
Notifications
You must be signed in to change notification settings - Fork 0
/
main.m
125 lines (88 loc) · 2.85 KB
/
main.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
%% simulate and estimate s(t)
data = csvread('elspot_prices.csv', 1, 1);
P_all = data(:,7); % DK1
N = 2*365;
P_all = P_all(1:2*N);
flength = length(P_all) - N; % forecast length
fplot = 731:761;
%% estimate ARIMA
arima_mdl = arima(1,1,1);
arima_fx = zeros(length(P_all)-N,3);
P = P_all(1:N);
[st, thetals] = estimate_st_ls(P);
X = P - st;
em = arima_mdl.estimate(X);
for i = N+1:length(P_all)
P = P_all(1:i-1);
st = simulate_st(thetals, i);
X = P - st(1:i-1);
[fx, ymse] = em.forecast(1, X);
fxn = fx - 1.96*sqrt(ymse);
fxp = fx + 1.96*sqrt(ymse);
arima_fx(i-N,1:3) = [fx fxn fxp] + st(end);
end
rmse_arima = sqrt(sum((P_all(N+1:end) - arima_fx(:,1)).^2)/flength);
figure('Name', 'ARIMA')
plot(731:761, P_all(731:761))
hold on
plot(731:761, arima_fx(1:31,1))
plot(731:761,arima_fx(1:31,2),'r:','LineWidth',2)
plot(731:761,arima_fx(1:31,3),'r:','LineWidth',2)
xlabel('Day')
ylabel('DKK/MWh')
legend('Real price', 'Forecasted price', '95% confidence')
%% estimate TV-ARFIMA
%lower = [0.05 -0.1 -0.2 -0.999 1];
%upper = [0.495 0.1 0.2 0.999 inf];
initial = [0.34 0 0.005 0.05 1];
lower = [0.05 -inf -inf -0.999 1];
upper = [0.495 inf inf 0.999 inf];
%initial = [0.25 0 0 0 1000];
%lower = [0.0001 -inf 0 -0.999 1];
%upper = [0.4999 inf inf 0.999 inf];
%initial = [0.4 0 0 1 1000];
tvarfima_fx = zeros(length(P_all)-N,3);
P = P_all(1:N);
[st, thetals] = estimate_st_ls(P);
X = P - st;
[d0,w,a,b,sigma_tv2,phi,theta,mu,sigma_arma2] = tvarfima_estimate(X, initial, lower, upper);
for i = N+1:length(P_all)
P = P_all(1:i-1);
st = simulate_st(thetals, i);
X = P - st(1:i-1);
[fx, fxn, fxp] = tvarfima_forecast(X, 1, d0,w,a,b,sigma_tv2,phi,theta,mu,sigma_arma2);
tvarfima_fx(i-N,1:3) = [fx fxn fxp] + st(end);
end
rmse_tvarfima = sqrt(sum((P_all(N+1:end) - tvarfima_fx(:,1)).^2)/flength)
%%
figure('Name', 'TV-ARFIMA')
plot(731:761, P_all(731:761))
hold on
plot(731:761, tvarfima_fx(1:31,1))
plot(731:761,tvarfima_fx(1:31,2),'r:','LineWidth',2)
plot(731:761,tvarfima_fx(1:31,3),'r:','LineWidth',2)
xlabel('Day')
ylabel('DKK/MWh')
legend('Real price', 'Forecasted price', '95% confidence')
%% estimate ARFIMA
arfima_fx = zeros(length(P_all)-N,1);
P = P_all(1:N);
[st, thetals] = estimate_st_ls(P);
X = P - st;
em = arfima_estimate(X, 'FWHI', [6 5]);
for i = N+1:length(P_all)
P = P_all(1:i-1);
st = simulate_st(thetals, i);
X = P - st(1:i-1);
fx = arfima_forecast(X, 1, em.d(1), em.AR, em.MA, em.mean, em.sigma2);
arfima_fx(i-N) = fx + st(end);
end
%rmse_arfima = sqrt(sum(((P_all(N+1:length(P_all)) - arfima_fx).^2)/(length(P_all)-N)));
plot(731:761, P_all(731:761))
hold on
plot(731:761, arfima_fx(1:31,1))
%plot(731:761,arfima_fx(1:31,2),'r:','LineWidth',2)
%plot(731:761,arfima_fx(1:31,3),'r:','LineWidth',2)
xlabel('Day')
ylabel('DKK/MWh')
legend('Real price', 'Forecasted price')