temp = load('temp.dat');

Temp = temp(:, 1);
Time = temp(:, 2);

co = cos(2*pi*Time/72);
si = sin(2*pi*Time/72);

[seasonb, seasonbint, seasonres] = regress(Temp, [ones(size(co)), co, si]);

x = zeros(size(seasonres, 1), 2);

lag = @(x, l) [NaN(l, 1); x(1:(end-l))];
for i = 1:2
    x(:, i) = lag(seasonres, i);
end

y = seasonres(3:end);
X = x(3:end, :);

[phi, phici, tempres] = regress(y, X);


subplot(2, 1, 1)
plot(Temp(3:end), '-k', 'LineWidth', 2)

set(gca, 'Box', 'off', 'FontSize', 20, ...
         'XLim', [0, 360], 'YLim', [36.6, 38.4], ...
         'XTick', 0:72:288, 'YTick', 36.5:0.5:38.5, ...
         'XTickLabel', 0:4, 'TickDir', 'out')
cS = strcat('Temperature (', char(176), 'C)');
xlabel('Time (days)', 'FontSize', 20)
ylabel(cS, 'FontSize', 20)

subplot(2, 1, 2)
plot(Temp(3:end), '-k', 'LineWidth', 2)
hold on;
plot(Temp(3:end)-seasonres(3:end)+tempres, '-r', 'LineWidth', 2)

set(gca, 'Box', 'off', 'FontSize', 20, ...
         'XLim', [0, 360], 'YLim', [36.6, 38.4], ...
         'XTick', 0:72:288, 'YTick', 36.5:0.5:38.5, ...
         'XTickLabel', 0:4, 'TickDir', 'out')
cS = strcat('Temperature (', char(176), 'C)');
xlabel('Time (days)', 'FontSize', 20)
ylabel(cS, 'FontSize', 20)

set(gcf, 'Position', [200, 100, 1400, 900])