% This is a script to generate ARMA polynomials and then simulated time series
% At the top of the script you can set NP=no. poles (same as no. zeros of \phi)
% and NZ=no. zeros (of \theta). Then you will use the mouse to move the cursor
% to where you want the poles to be. This program assumes that if the polynomial is
% of even order all poles (or zeros) will occur in complex conjugate pairs. If NP or
% NZ are odd, it assumes the even part will occur in conjugate pairs. If you choose
% NP=1 and NZ=0, the pole must be on the real axis so it will place the pole at the
% x component of thepoint you select. If you set NP=2 and NZ=1, you will first be 
% asked to input the pole locations. This requires you to select only one point for
% two poles because the other one is set to be the complex conjugate of the one you chose.
% Next you will be asked to input a zero and since the order is  one, it again takes the 
% real part. If you set NP=3, for example, you first will be asked to place the single
% pole on the real axis.

NP=2;		% no. poles (zeros of phi)
NZ=0;		% no. zeros (zeros of theta)
NSAMP=500;  % no. samples to generate
INIT=0;     % no samples to discard before plotting
clear phi theta;
[phi,theta]=pickpoly(NP,NZ,pi/6.0,0); % in these polys the constant term is 1st
% NOTE the phi,theta produced are in matlab format, IE coeffs of highest power first
phi=phi/phi(1); % normalize to make first coeff 1
theta=theta/theta(1);
phiAR=fliplr(phi);phiAR=phiAR/phiAR(1); % make leading coeff 1
thetaMA=fliplr(theta);thetaMA=thetaMA/thetaMA(1);
disp(['\phi coefficients in AR form: ',num2str(phiAR)]);
disp(['\theta coefficients in MA form: ',num2str(thetaMA)]);

%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%% plot the 1st NSAMP values of the impulse response
figure;
x=zeros(1,NSAMP);
x(1)=1;				% make an impulse
y=filter(theta,phi,x); % filter(B,A,X) is AY=BX so A is AR and B is MA
subplot(211);
plot(y);
title('impulse response : shows the filter weights');
%%%%%%%%%% plot the NSAMP values of an ARMA sequence
subplot(212);
n=NSAMP+INIT;
x=randn(n,1);
y=filter(theta,phi,x); % filter(B,A,X) is AY=BX so A is AR and B is MA
out=y(INIT+1:INIT+NSAMP);
plot(y(INIT+1:INIT+NSAMP));
v=axis;
text(v(1)+.2*(v(2)-v(1)),v(3)+.8*(v(4)-v(3)),['theta:',...
      sprintf(' %4.2f',theta),sprintf('\n phi:'),sprintf(' %4.2f',phi)]);
title([int2str(NSAMP),' samples of ARMA(',int2str(NP),',',int2str(NZ),') simulation']);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%  plot |W(\lambda)|^2
NLAM=200;	% no. pts to plot in  [0,2\pi]
dlam=2*pi/NLAM;
lams=[0:dlam:2*pi-dlam];
e=exp(-i*lams);

thetaspec=polyval(theta,e);
phispec=polyval(phi,e);

W=thetaspec./phispec;

%%%%%%%%%%%%%%  first using linear scale
figure;
plot(lams,(abs(W).^2)/(2*pi));
xlabel(' 0 \leq \lambda < 2\pi');
v=axis;
text(v(1)+.2*(v(2)-v(1)),v(3)+.8*(v(4)-v(3)),['theta:',...
      sprintf(' %4.2f',theta),sprintf('\n phi:'),sprintf(' %4.2f',phi)]);
title('spectral density - linear scale');

%%%%%%%%%%%%%%  next using log scale
hlogp=figure;
semilogy(lams,(abs(W).^2)/(2*pi));
xlabel(' 0 \leq \lambda < 2\pi');
v=axis;
text(v(1)+.2*(v(2)-v(1)),v(3)+.8*(v(4)-v(3)),['theta:',...
      sprintf(' %4.2f',theta),sprintf('\n phi:'),sprintf(' %4.2f',phi)]);
title('spectral density - log scale');

%%%%%%%%%%%%%%  theoretical spectral DF
figure;
plot(lams,cumsum(abs(W).^2));	% cumsum is like an integral;
xlabel(' 0 \leq \lambda < 2\pi');
v=axis;
text(v(1)+.2*(v(2)-v(1)),v(3)+.8*(v(4)-v(3)),['theta:',...
      sprintf(' %4.2f',theta),sprintf('\n phi:'),sprintf(' %4.2f',phi)]);
title('spectral CDF - linear scale');



%%%%%%%%%%%%%% next periodogram 

WINLEN=16;
win=window(@hanning,WINLEN);    % use hanning type window; others possible, type help window
win=win/(sum(win));             % normalize to have integral (sum) of 1

pgx=abs(fft(out));
n=length(out);
pgx=(pgx.*pgx)/(2*pi*n);
if WINLEN > 0
    pgx=conv(pgx,win);                      % convolves pgx with win
    pgx=pgx(WINLEN:length(pgx) -WINLEN);    % cuts off the ends
end
figure;
semilogy(pgx(2:end));
title(['smoothed periodogram of ARMA(',int2str(NP),',',int2str(NZ),') simulation']);
 

gammax=acf(out,0);
figure;
plot(gammax(1:floor(n/4)));
title(['ACF of ARMA(',int2str(NP),',',int2str(NZ),') simulation']);