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

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

Answer 1

[Mathieu NOE]

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');

[Ikmal Rosli]

Thanks Mathieu!

[Mathieu NOE]

You're welcome