#!/usr/bin/env perl
#
#
# ssa.pl filename [maxlag] [arxord] [predhorizon]
#
# Study a system y=f(x) using an ARX model
#
# filename points to a text file in which each line contains
# two numbers separated by white space.  The first line contains 
# x[t]  y[t] the next x[t+1] y[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
#
# prediction horizon is 0 by default
#
# The output is a color eps file "filename.fig.eps" containing a ts plot,
# the acf plot, a transfer function estimate 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 "ssa.pl filename [maxlag] [arxorder] [predhorizon]\n";
    exit(-1);
}

$filename = $ARGV[0];
if ($#ARGV>0) { 
    $maxlag = $ARGV[1];
} else {
    $maxlag = 100;
}
if ($#ARGV>1) { 
    $arxord = $ARGV[2];
} else {
    $arxord = 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(:,1);\n";
print MATIN "y=matlabinput(:,2);\n";
print MATIN "z(:,1)=x; z(:,2)=y;\n";
print MATIN "n = length(x);\n"; 
print MATIN "subplot(2,2,1);\n";
print MATIN "plot3(1:n,x,y);\n";
#print MATIN "set(gca,'FontSize',14);\n";
print MATIN "xlabel('time');\n";
print MATIN "ylabel('input');\n";
print MATIN "zlabel('output');\n";
print MATIN "t = sprintf('input/output time series of $filename');\n";
print MATIN "title(t);\n";
print MATIN "acfall=xcov(x,y,$maxlag,'coeff');\n";
print MATIN "acf = acfall($maxlag+1:2*$maxlag+1);\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('XCORR(input,output)');\n";
print MATIN "t = sprintf('XCORR of $filename (%g %% sig at p=0.05)',sigp);\n";
print MATIN "title(t);\n";
print MATIN "s = spa(z);\n";
print MATIN "subplot(2,2,4);\n";
print MATIN "loglog(s(2:length(s),2));\n";
print MATIN "xlabel('log(freq)');\n";
print MATIN "ylabel('log(abs(xfer))');\n";
print MATIN "t = sprintf('transfer function amplitude  of $filename');\n";
print MATIN "title(t);\n";
print MATIN "th=arx(z(1:n/2,:),[$arxord 1 0]);\n";
print MATIN "yp=predict(z(n/2+1:n,:),th,$predhorizon);\n";
print MATIN "subplot(2,2,3);\n";
print MATIN "plot3(1:n,x,y,n/2+1+$arxord:n,x(n/2+1+$arxord:n),yp($arxord+1:length(yp)),'r--');\n";
print MATIN "xlabel('time');\n";
print MATIN "ylabel('input');\n";
print MATIN "zlabel('output/prediction');\n";
print MATIN "t = sprintf('$predhorizon step ahead i/o preds of ARX($arxord) 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 xcfs 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";