Hi Frank,
The code seems correct as you have correctly implemented the theoretical autocovariance function of an AR(1) process and used it to calculate the theoretical spectral function. You have also correctly implemented the Fast Fourier Transform to calculate the spectral dsitance of time series.
If you want to interpret the data as quaterly and want the units of the x-axis as "Cycles per year" or "Periods per Cycle", you would need to adjust the frequency accordingly.
For "cycles per year", you can divide the frequencies by the 4. Then we can adjust the x-axis values to match the quaterly data. Here is the code for plotting this:
freq_cycles_per_year = (0:N-1) / (T/4);
x_adj = x/4; % Adjust x for quarterly data
figure
plot(x_adj, spec_theo(1:N))
hold on
plot(x_adj, spec_fft(2:N+1)) %cutoff the first value
legend('theoretical','fft')
Similarly, for "periods per cycle", we can take the reciprocal of the frequencies to get periods. Here is the code for the same:
periods_per_cycle = 1 ./ freq_cycles_per_year;
periods_per_cycle(1) = NaN; % Avoid division by zero for the first element
x_adj = 4./x; % Adjust x for quarterly data
figure
plot(x_adj,spec_theo(1:N))
hold on
plot(x_adj,spec_fft(2:N+1)) %cutoff the first value
legend('theoretical','fft')
Hope this help!
