clear;
%pathSSH='d:\TSP\LR4'
pathSSH='~/2024/lectures/TSP/LR4'

N_signal=1024;

% ARMA process generating
arnoise(1)=0.5*randn(1);
arnoise(2)=0.3*arnoise(1)+0.5*randn(1);

for (i=3:1:N_signal)
    arnoise(i)=0.3*arnoise(i-1)+0.7*arnoise(i-2)+0.5*randn(1);
end;
Noisestd=1;
arnoise=Noisestd*arnoise;

plot(arnoise)
ar_order=3;
 [arm,ref1] = ar(arnoise,ar_order,'burg','now') ;
arm.a(1)
a1=-arm.a(2);
a2=-arm.a(3);
nsko=sqrt(arm.NoiseVariance);


% ARMA process generating
armodel(1)=nsko*randn(1);
armodel(2)=a1*arnoise(1)+nsko*randn(1);

for (i=3:1:N_signal)
    armodel(i)=a1*armodel(i-1)+a2*armodel(i-2)+nsko*randn(1);
end;
Noisestd=1;
armodel=Noisestd*armodel;

plot(1:N_signal,arnoise,1:N_signal,armodel)

% file reading
filename_SSH=sprintf('%s/lambda033_phi34pc__all.dat',pathSSH);

 fin=fopen(filename_SSH,'rt');
 Adat=fscanf(fin,'%f',[2 inf]);
 fclose(fin);

 dates=Adat(1,:);
 x=Adat(2,:);
 l=size(Adat);

 N=l(2);
 
 
plot(dates,x);

 signal_centered=x-mean(x);
 
for(tau=1:1:N)
 acf(tau)=0;
 for(j=1:1:N-tau)
    acf(tau)=acf(tau)+signal_centered(j)*signal_centered(j+tau-1);
 end;
  acf(tau)= acf(tau)/(N);
end;
plot(acf);



spectr_dens=fft(acf);
plot(abs((spectr_dens)))



%array of dates to predict
 m=mean(x);
 dt=dates(2)-dates(1);
 dates_to_predict=dates(N):dt:dates(N)+0.8
 
 % remove polynomial trend
poly=1;
deg=3;
if(poly)
 [x_without_trend,poly_pred]=predict_poly(dates,x,dates_to_predict, deg);
else
    x_without_trend=x;
end;
 
plot(dates,x,dates,x-x_without_trend,dates_to_predict,poly_pred)


 [ spectr, freq] = spect_fftn(dates, x_without_trend)
 plot(freq , abs(spectr))

%harmonic trend
  Periods=[1 0.5 14];
  [Xsinh, harm_pred] = predict_harm(dates,x_without_trend,Periods,dates_to_predict);
    plot(dates,x,dates_to_predict,harm_pred+poly_pred,'red',dates,x-Xsinh+m);  
plot(dates,Xsinh)
     
  %autoregression
   ar_order=25;
  [XAR_pred,cf] = predict_ar(dates,Xsinh,dates_to_predict,ar_order);
  plot(dates, Xsinh,dates_to_predict,XAR_pred)
 
  
 %add all prediction together
 dates2=dates_to_predict;
 z=harm_pred+poly_pred+XAR_pred;
 l=size(z);

 N2=l(2);
 
 %plot and output to file
 plot(dates,x,dates_to_predict,z);  

  foutname= sprintf('%s/SSH_predition.dat',pathSSH);
  fout=fopen(foutname, 'wt');
  for (j=1:1:N2)
    fprintf(fout,'%6.2f %10.8e\n',dates_to_predict(j),z(j));
  end; 
  fclose(fout);
  