I'm calculating the harmonics of a time series. How do I determine the significance of the harmonic fit onto the data?
12 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Ikmal Rosli
le 5 Juil 2021
Commenté : Mathieu NOE
le 6 Juil 2021
I am calculating the the first few harmonics of a time series. The data and code below is a replication of an example from a textbook from Wilks (2011), pg 371 - 381, Statistical Methods in the Atmospheric Sciences. I've got the amplitude and phase shift right. But I'm not sure how to evaluate the significance of the fit.
clear;clc
% fitting a harmonic
yt = [22.2 22.7 32.2 44.4 54.8 64.3 68.8 67.1 60.2 49.5 39.3 27.4]'; % data
t = (1:12)';
n = 12;
% first harmonic
omegaA1 = cosd(360*t/n);
omegaB1 = sind(360*t/n);
% second harmonic
omegaA2 = cosd(360*t*2/n);
omegaB2 = sind(360*t*2/n);
% regress
b1 = [ones(length(yt), 1) omegaA1 omegaB1 omegaA2 omegaB2]\yt;
% get first harmonic and phase shift;
A1 = hypot(b1(2), b1(3)); % first harmonic amplitude
if b1(2) > 0
phi1 = atand(b1(3)/b1(2))
elseif b1(2) < 0
phi1 = atand(b1(3)/b1(2));
if phi1 < 180
phi1 = phi1 + 180;
elseif phi1 > 180
phi1 = phi1 -180;
end
elseif b(1) == 0
phi1 = 90;
end
% second harmonic and phase shift;
A2 = hypot(b1(4), b1(5)); % second harmonic amplitude
if b1(4) > 0
phi2 = atand(b1(5)/b1(4));
elseif b1(4) < 0
phi2 = atand(b1(5)/b1(4));
if phi2 < 180
phi2 = phi2 + 180;
elseif phi2 > 180
phi2 = phi2 - 180;
end
elseif b1(4) == 0
phi2 = 90
end
% fit and plot
close all
yfit1 = b1(1) + A1*cosd(360*t/n - phi1);
plot(t, yt, 'o')
hold on
plot(t, yfit1)
hold off
[~, p] = corr(yt,yfit1); % I am not quite sure about this part
0 commentaires
Réponse acceptée
Mathieu NOE
le 5 Juil 2021
hello
why not simply compute the mean squared error between your fitted data and the input data
I also tried including the second harmonic , but it did not improve the fit...
clearvars;clc
% fitting a harmonic
yt = [22.2 22.7 32.2 44.4 54.8 64.3 68.8 67.1 60.2 49.5 39.3 27.4]'; % data
t = (1:12)';
n = 12;
% first harmonic
omegaA1 = cosd(360*t/n);
omegaB1 = sind(360*t/n);
% second harmonic
omegaA2 = cosd(360*t*2/n);
omegaB2 = sind(360*t*2/n);
% regress
b1 = [ones(length(yt), 1) omegaA1 omegaB1 omegaA2 omegaB2]\yt;
% get first harmonic and phase shift;
A1 = hypot(b1(2), b1(3)); % first harmonic amplitude
if b1(2) > 0
phi1 = atand(b1(3)/b1(2))
elseif b1(2) < 0
phi1 = atand(b1(3)/b1(2));
if phi1 < 180
phi1 = phi1 + 180;
elseif phi1 > 180
phi1 = phi1 -180;
end
elseif b(1) == 0
phi1 = 90;
end
% second harmonic and phase shift;
A2 = hypot(b1(4), b1(5)); % second harmonic amplitude
if b1(4) > 0
phi2 = atand(b1(5)/b1(4));
elseif b1(4) < 0
phi2 = atand(b1(5)/b1(4));
if phi2 < 180
phi2 = phi2 + 180;
elseif phi2 > 180
phi2 = phi2 - 180;
end
elseif b1(4) == 0
phi2 = 90
end
% fit and plot
close all
yfit1 = b1(1) + A1*cosd(360*t/n - phi1);
yfit2 = b1(1) + A1*cosd(360*t/n - phi1)+ A2*cosd(360*t/n - phi2);
mse1 = sqrt(mean((yfit1-yt).^2)) % mean squared error for fit #1
mse2 = sqrt(mean((yfit2-yt).^2)) % mean squared error for fit #2
plot(t, yt, 'o',t, yfit1,t, yfit2)
legend('data','fit 1' ,'fit 2');
2 commentaires
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur Hypothesis Tests dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!