#!/usr/bin/env perl
#
#
# tsa.pl filename [maxlag] [arord] [predictionhorizon]
#
# Study a time series x[t] using AR models
#
# filename points to a text file in which each line contains
# a single number.  The first line contains x[t], the next x[t+1], and
# so on.
# 
# maxlag is the maximum lag for the autocorrelation the default maxlag is 100
#
# arord is the order of ar model fit to the first 1/2 of the sequence and
# used to predict the second half.  default is 32
#
# predictionhorizon is for how far ahead prediction are made - default 1
#
# The output is a color eps file "filename.fig.eps" containing a ts plot,
# the acf plot, a loglog fft plot, and a plot of the sequence with the 
# predictions for the second half based on the model fit on the first half
#
#
# The program will also state how many of the >0 lag components 
# are significant
#
#

use IPC::Open2;
use FileHandle;

if ($#ARGV<0 || $#ARGV>3) {
    print STDERR "tsa.pl filename [maxlag] [arorder] [predhorizon]\n";
    exit(-1);
}

$filename = $ARGV[0];
if ($#ARGV>0) { 
    $maxlag = $ARGV[1];
} else {
    $maxlag = 100;
}
if ($#ARGV>1) { 
    $arord = $ARGV[2];
} else {
    $arord = 32;
}
if ($#ARGV>2) { 
    $predhorizon = $ARGV[3];
} else {
    $predhorizon = 1;
}


system "cp $filename matlabinput.txt";

open2(MATOUT,MATIN,"matlab -display none");
MATOUT->autoflush(1);
MATIN->autoflush(1);

print MATIN "load('matlabinput.txt');\n";
print MATIN "x=matlabinput;\n";
print MATIN "acfall=xcov(x,$maxlag,'coeff');\n";
print MATIN "acf = acfall($maxlag+1:2*$maxlag+1);\n";
print MATIN "n = length(x);\n"; 
print MATIN "bound95 = 1.96./sqrt(n:-1:n-$maxlag);\n";
print MATIN "lags = 0:$maxlag;\n";
print MATIN "sigp = 100*sum(abs(acf(2:$maxlag))>bound95(2:$maxlag)')/$maxlag\n";
print MATIN "subplot(2,2,2);\n";
print MATIN "plot(lags,acf,lags,bound95,'r--',lags,-bound95,'r--');\n";
#print MATIN "set(gca,'FontSize',14);\n";
print MATIN "xlabel('Lag');\n";
print MATIN "ylabel('Acf');\n";
print MATIN "t = sprintf('ACF of $filename (%g %% sig at p=0.05)',sigp);\n";
print MATIN "title(t);\n";
print MATIN "subplot(2,2,1);\n";
print MATIN "plot(x);\n";
#print MATIN "set(gca,'FontSize',14);\n";
print MATIN "xlabel('time');\n";
print MATIN "ylabel('signal');\n";
print MATIN "t = sprintf('time series of $filename');\n";
print MATIN "title(t);\n";
print MATIN "f = abs(fft(x));\n";
print MATIN "f = f(1:n/2);\n";
print MATIN "subplot(2,2,4);\n";
print MATIN "loglog(f);\n";
print MATIN "xlabel('log(freq)');\n";
print MATIN "ylabel('log(abs(fft))');\n";
print MATIN "t = sprintf('PSD of $filename');\n";
print MATIN "title(t);\n";
print MATIN "th=ar(x(1:n/2),$arord);\n";
print MATIN "yp=predict(x(n/2+1:n),th,$predhorizon);\n";
print MATIN "subplot(2,2,3);\n";
print MATIN "plot(1:n,x,n/2+1+$arord:n,yp($arord+1:length(yp)),'r--');\n";
print MATIN "xlabel('time');\n";
print MATIN "ylabel('signal/prediction');\n";
print MATIN "t = sprintf('$predhorizon step ahead preds of AR($arord) on $filename');\n";      
print MATIN "title(t);\n";
print MATIN "print -depsc $filename.figs.eps\n";
print MATIN "quit;\n";

@linesout = <MATOUT>;
close(MATOUT);
close(MATIN);

$globout = join("",@linesout);
#print $globout;
$globout =~ /sigp =.*\n.*\n(.*)\n/;
$sigp = $1;
$sigp =~ s/\s+//g;

print "Figures written to $filename.figs.eps\n";
print "$sigp % of the acfs are significant at p=0.05 (expect less than 5 %)\n";
if ($sigp>5) {
    print "This sequence may have some (linear) predictability!\n";
} else {
    print "This sequence is probably not linearly predictable\n";
}

system "rm -f matlabinput.txt";