Spectral Analysis with fft

Hi guys,

currently, I am trying to understand how spectral analysis of time series can be done in matlab, but unfortunately I have some difficulties in understanding the concept.

As a first stept I would like to compare the theoretical spectral density of an AR(1)-process with the one implied by a realization of an AR(1)-process.

This is what I have so far

rng(786786)
AR1 = arima('Constant',0,'AR',{0.9},'Variance',0.01^2);
Y1 = simulate(AR1,1000);

ts=Y1(end-200+1:end); %Take last 200 observations
T = length(ts); %Length of time series

%Calculate theoretical spectral function 
k=1:1:T-1;
theo_gamma=[0.01^2/(1-0.9^2) (0.9.^k)*0.01^2/(1-0.9^2)]'; %theoretical autocovariance fuction of an AR(1)-process
omega=2*pi*k'./T;
spec_theo=zeros(1,length(k));

for j=k
    spec_theo(j)=theo_gamma(1)/(2*pi)+1/pi*theo_gamma(2:end)'*cos(omega(j)*k');
end

%Calculating spectral density with fft
fft_ts = fft(ts);
spec_fft = abs(fft_ts).^2/(2*pi*T);

N=floor(T/2)+1;
x=1:N;

figure
plot(x,spec_theo(1:N))
hold on
plot(x,spec_fft(2:N+1)) %cutoff the first value
legend('theoretical','fft')

I divided the the period length of 2*pi into equidistant subintervals according to the time series length. My first question: Is the code correct?

Further, I would like to analyze some real data, so imagine that the 200 Realizations of the AR(1)-process correspond to 50 Years of quaterly data. In some textbooks, the are graphs showing the units of the x-Axis as "Cycles per year" or "Periods per Cycle". What do I have to do in order to get such graphs in my example, if we interpret the data as quaterly?

Thanks a lot!

Frank

Answer 1

[Jaynik]

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!