I'm calculating the harmonics of a time series. How do I determine the significance of the harmonic fit onto the data?

16 views (last 30 days)
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

Accepted Answer

Mathieu NOE
Mathieu NOE on 5 Jul 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');

More Answers (0)

Products


Release

R2021a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!