getFred = function(file) { cat("reading", file, "...");
        skip = 0;
        numEmpty = 0;
        repeat {
                a = scan(file, skip = skip, nlines = 1, what = "",
                         sep = "\n");
                if (length(a) == 0) {
                        if (numEmpty > 0)
                                stop("Didn't find the DATE line");
                        numEmpty = numEmpty + 1;
                } else {
                        if (substr(a, 1, 4) == "DATE")
                                break;
                }
                skip = skip + 1;
        }
        a = read.table(file, skip = skip, header = TRUE);
        cat("done\n");
        if (!("VALUE" %in% names(a)))
                stop(paste("Bad skip in file", file, ":", skip));
        b = substring(a[1, 1], 1, 4);
        c = substring(a[1, 1], 6, 7);
        ts(a[,2], frequency = 12, start = as.numeric(c(b, c)));
};

cPlot = function(x, text = names(x), l) {
        if (missing(l))
                l = max(Mod(x));
        l = c(-1, 1) * l;
        plot(c(0, 1.1 * x), type = "n", asp = 1, xlim = l, ylim = l);
        arrows(0, 0, Re(x), Im(x));
        text(1.05 * x, text);
};
xyLim = 0.1;

fredPlot = function(x, l, f = fred)
                plot(ts(x, frequency = frequency(f), start = start(f)),
                     ylab = l);

BUSINV   = getFred("BUSINV.txt");
CPIAUCSL = getFred("CPIAUCSL.txt");
CPILFESL = getFred("CPILFESL.txt");
GS10     = getFred("GS10.txt");
GS1      = getFred("GS1.txt");
HOUST    = getFred("HOUST.txt");
INDPRO   = getFred("INDPRO.txt");
PAYEMS   = getFred("PAYEMS.txt");
PPIFGS   = getFred("PPIFGS.txt");
PPILFE   = getFred("PPILFE.txt");
TCU      = getFred("TCU.txt");
UNRATE   = getFred("UNRATE.txt");

fred = ts.intersect(
                diff(log(BUSINV)),
                diff(log(CPIAUCSL)),
                diff(log(CPILFESL)),
                GS10,
                GS1,
                diff(log(HOUST)),
                diff(log(INDPRO)),
                diff(log(PAYEMS)),
                diff(log(PPIFGS)),
                diff(log(PPILFE)),
                TCU,
                UNRATE
        );
shortNames = c("BU", "CA", "CL", "1", "10", "HO", "IN", "PY", "PA",
               "PL", "TC", "UN");

fred = ts(apply(fred, 2, function(x) (x - mean(x)) / sd(x)),
          frequency = frequency(fred), start = start(fred));

# default 10% tapering:
fredT = spec.taper(fred);

n = nrow(fredT);
fredF = apply(fredT, 2, fft);
f = (0:(n - 1)) / n;
flo = 1/120;
fhi = 1/12;
fredF[f <= flo | f >= fhi,] = 0;
fredC = apply(fredF, 2, fft, inv = TRUE) / n;

fredS = svd(fredC);

graphics.off();
cPlot(fredS$v[,1], dimnames(fred)[[2]]);

x11();
par(mfrow = c(2, 2));

fredPlot(Re(fredS$u[,1]), "Real");
fredPlot(Im(fredS$u[,1]), "Imag");
fredPlot(Mod(fredS$u[,1]), "Mod");
fredPlot(unwrap(Arg(fredS$u[,1]) / (2 * pi)), "Arg");

if (!exists("plotCsvd"))
        plotCsvd = TRUE;
if (!exists("plotPause"))
        plotPause = FALSE;
if (plotCsvd) {
        if (!plotPause)
                scan();
        x11();
        for (i in 1:nrow(fred)) {
                if (plotPause)
                        scan();
                cPlot(fredS$u[i, 1] * Conj(fredS$v[,1]), shortNames, xyLim);
                a = time(fred)[i];
                b = floor(a);
                c = 12 * (a - b);
                title(paste(month.abb[1 + round(c)], b));
        }
        scan();
}

# look at 2nd component
x11();
cPlot(fredS$v[,2], dimnames(fred)[[2]]);

x11();
par(mfrow = c(2, 2));

fredPlot(Re(fredS$u[,2]), "Real");
fredPlot(Im(fredS$u[,2]), "Imag");
fredPlot(Mod(fredS$u[,2]), "Mod");
fredPlot(unwrap(Arg(fredS$u[,2]) / (2 * pi)), "Arg");

if (!exists("plotCsvd2"))
        plotCsvd2 = TRUE;
if (plotCsvd2) {
        if (!plotPause)
                scan();
        x11();
        for (i in 1:nrow(fred)) {
                if (plotPause)
                        scan();
                z = fredS$d[1] * fredS$u[i, 1] * Conj(fredS$v[,1]) + 
                    fredS$d[2] * fredS$u[i, 2] * Conj(fredS$v[,2]);
                cPlot(z, shortNames, 1);
                a = time(fred)[i];
                b = floor(a);
                c = 12 * (a - b);
                title(paste(month.abb[1 + round(c)], b));
        }
        scan();
}

# approximate the real and imaginary parts of fredS$u[i, 1] with fredT:
fredL = lm(cbind(Re(fredS$u[, 1]), Im(fredS$u[, 1])) ~ ., fredT);
summary(fredL);
u1 = ts(predict(fredL, fred),
        start = start(fred),
        frequency = frequency(fred));
u1c = complex(real = u1[, 1], imaginary = u1[, 2]);

x11();
par(mfrow = c(2, 2));
fredPlot(Re(u1c), "Real");
fredPlot(Im(u1c), "Imag");
fredPlot(Mod(u1c), "Mod");
fredPlot(unwrap(Arg(u1c) / (2 * pi)), "Arg");

if (!exists("plotPredict"))
        plotPredict = TRUE;
if (plotPredict) {
        if (!plotPause)
                scan();
        x11();
        for (i in 1:nrow(fred)) {
                if (plotPause)
                        scan();
                cPlot(u1c[i] * Conj(fredS$v[,1]), shortNames, xyLim);
                a = time(fred)[i];
                b = floor(a);
                c = 12 * (a - b);
                title(paste(month.abb[1 + round(c)], b));
        }
        scan();
}

u1l = complex(real      = lowess(u1[, 1], f = 0.06)$y,
              imaginary = lowess(u1[, 2], f = 0.06)$y);
u1l = ts(u1l,
        start = start(fred),
        frequency = frequency(fred));

x11();
par(mfrow = c(2, 2));
fredPlot(Re(u1l), "Real");
fredPlot(Im(u1l), "Imag");
fredPlot(Mod(u1l), "Mod");
fredPlot(unwrap(Arg(u1l) / (2 * pi)), "Arg");

if (!exists("plotLowess"))
        plotLowess = TRUE;
if (plotLowess) {
        if (!plotPause)
                scan();
        x11();
        for (i in 1:nrow(fred)) {
                if (plotPause)
                        scan();
                cPlot(u1l[i] * Conj(fredS$v[,1]), shortNames, xyLim);
                a = time(fred)[i];
                b = floor(a);
                c = 12 * (a - b);
                title(paste(month.abb[1 + round(c)], b));
        }
        scan();
}