temp = load('temp.dat'); Temp = temp(:, 1); Time = temp(:, 2); co = cos(2*pi*Time/72); si = sin(2*pi*Time/72); [tempb, tempbint, tempres] = regress(Temp, [ones(size(co)), co, si]); orderP = 22; n = 330; x = zeros(size(tempres, 1), orderP); lag = @(x, l) [NaN(l, 1); x(1:(end-l))]; for i = 1:22 x(:, i) = lag(tempres, i); end y = tempres((orderP+1):end); X = x((orderP+1):end, :); phi = regress(y, X); tacf = conv(flipud(tempres), tempres); acf = tacf(size(tempres, 1):(size(tempres, 1)+orderP))./tacf(size(tempres, 1)); pacf = zeros(orderP, 1); pacf(1) = acf(2); for i = 2:orderP acfrec = toeplitz(acf(1:i)); acfcur = acf(2:(i+1)); tpacf = acfrec\acfcur; pacf(i) = tpacf(i); end subplot(3, 1, 1) plot(phi, 'ok', 'MarkerSize', 6) hold on; plot(phi, ':k') plot([0, orderP+1], [0, 0], '-k') set(gca, 'FontSize', 20, 'XLim', [-0.5, 23.2], 'XTick', 0:5:20, 'YLim', [-0.4, 1.1], 'YTick', -0.5:0.5:1, 'box', 'off') xlabel('Lag', 'FontSize', 20) ylabel('$\phi_i$', 'FontSize', 20, 'Interpreter', 'LaTex') tx = [0:orderP; 0:orderP]; ty = [zeros(1, orderP+1); acf']; subplot(3, 1, 2) plot(tx, ty, '-k', 'LineWidth', 2) hold on; plot([0, orderP+1], [0, 0], '-k') plot([0, orderP+1], [1.96/sqrt(size(tempres, 1)), 1.96/sqrt(size(tempres, 1))], ':k') plot([0, orderP+1], [-1.96/sqrt(size(tempres, 1)), -1.96/sqrt(size(tempres, 1))], ':k') set(gca, 'FontSize', 20, 'XLim', [-0.5, 23.2], 'XTick', 0:5:20, 'YLim', [-0.2, 1.1], 'YTick', -0.5:0.5:1, 'box', 'off') xlabel('Lag', 'FontSize', 20) ylabel('ACF', 'FontSize', 20) tx = [1:orderP; 1:orderP]; ty = [zeros(1, orderP); pacf']; subplot(3, 1, 3) plot(tx, ty, '-k', 'LineWidth', 2) hold on; plot([0, orderP+1], [0, 0], '-k') plot([0, orderP+1], [1.96/sqrt(size(tempres, 1)), 1.96/sqrt(size(tempres, 1))], ':k') plot([0, orderP+1], [-1.96/sqrt(size(tempres, 1)), -1.96/sqrt(size(tempres, 1))], ':k') set(gca, 'FontSize', 20, 'XLim', [-0.5, 23.2], 'XTick', 0:5:20, 'YLim', [-0.5, 1], 'YTick', -0.5:0.5:1, 'box', 'off') xlabel('Lag', 'FontSize', 20) ylabel('pACF', 'FontSize', 20) set(gcf, 'Position', [200, 0, 1400, 1000])