####################################################################
## Purpose:
## Predict SBs due to phasing out of SNAs for LPI.
##
## Execute following on dos command prompt (omit the "##"):
## "S:\\R\\rw2001\\bin\\Rterm.exe" -q --no-restore --no-save < LPI_SBs.r 
## Unix:
## /usr/bin/R -q --no-restore --no-save < LPI_SBs.r 
####################################################################

options(echo = FALSE)

#---------------------------------------------------------------------
# SN mvt data in %, read/store original data and compute mvts,
# also compute NSN mvts:
#           3/03   4/03  1/04  2/04 3/04  4/04  1/05  2/05  3/05  4/05  1/06 
#            1      2     3     4    5     6     7     8     9    10    11
xlabels <- c("Sep/'03", "Mar/'04", "Sep/'04", "Mar/'05", "Sep/'05", "Mar/'06");

#ALL:
#MvtSN <- c(0.155, 0.075, 0.009, 0.004, 0.145, 0.094, 0.013, 0.004, 0.13, 0.059, 0.004);
#PRIV:
#MvtSN <- c(0.203, 0.096, 0.011, 0.005, 0.191, 0.117, 0.016, 0.006, 0.175, 0.078, 0.005);
#PUB:
MvtSN <- c(0.007, 0.010, 0.001, 0.000, 0.005, 0.026, 0.005, 0.000, 0.003, 0.006, 0.001);

x1 <- scan("PUBINDUST_orig.txt", list(0,0,0));
xd1 <- x1[[3]];

startorigindex <- 24;
MvtTot <- 100*diff(xd1[startorigindex:length(xd1)],1)/
          xd1[startorigindex:(length(xd1)-1)];
MvtNSN <- MvtTot - MvtSN;

# read/store corresponding X11 trend derived using original to plot below:
t1 <- scan("PUBINDUST_trend.txt", list(0,0,0));
talldata1 <- log(t1[[3]]);
td1 <- talldata1[startorigindex:length(talldata1)];

#---------------------------------------------------------------------
# Plot all movements:

#bitmap(file="Mvts_allindustries.png", width=5, height=5, res=400)
png(file="Mvts_allindustries.png", width = 680, height = 450, pointsize = 12, bg="white", res = 200);

startyear <- 2003;
startperiod <- 3;
freq <- 4;
ts1 <- ts( cbind(MvtTot, MvtSN, MvtNSN), start=c(startyear, startperiod), frequency=freq );
ts.plot(ts1, gpars=list(ylab="% Change from previous period", col=c(1,2,4), frame = TRUE,
        axes=FALSE), main="% mvts in original for Safety Net Adjustments");
ts1;

legend(2005.5, 0.6,
       c("Total","SN jobs","Not-SN jobs"),
       col=c("black","red","blue"),
       text.col="black",
       lty=c(1,1,1),
       lwd=c(1,1,1),
       pch=c(" ", " ", " "),
       merge=FALSE,
       bty='n');

axis(1, axTicks(1), labels=xlabels, las=1)
axis(2);

#---------------------------------------------------------------------
# compute log of original; reconstruct log(O_nsn) vs. t:

o <- log(xd1[startorigindex:(length(xd1))]);

o_nsn <- 0;

# starting value in NSN series is arbitrary in order to reconstruct o_nsn.
# The residual after trend fitting is what matters for seas. factor estimation:

o_nsn[[1]] <- 4.575; #o[[1]];
for( t in 2:length(o) ) {
   o_nsn[[t]] <- o_nsn[[1]] + sum( log(1 + (MvtNSN[1:(t-1)]/100)) );
}

#---------------------------------------------------------------------
# fit trends to log of originals:

t_year <- 2003.0 + 0.25*cumsum(rep(1,length(o)));
trend_o <- lm(o ~ t_year + 1);
trend_o_nsn <- lm(o_nsn ~ t_year + 1);

par_o <- coefficients(summary(trend_o));
inter_o <- par_o[1,1];
slope_o <- par_o[2,1];
resid_o <- residuals(trend_o);

par_o_nsn <- coefficients(summary(trend_o_nsn));
inter_o_nsn <- par_o_nsn[1,1];
slope_o_nsn <- par_o_nsn[2,1];
resid_o_nsn <- residuals(trend_o_nsn);

#---------------------------------------------------------------------
# plot log trends and original data.

#bitmap(file="LogTrends_and_originals.png", width=5, height=5, res=400)
png(file="LogTrends_and_originals.png", width = 680, height = 450, pointsize = 12, bg="white", res = 200);

startyear <- 2003;
startperiod <- 2;
freq <- 4;
ts3 <- ts( cbind(o,o_nsn,td1), start=c(startyear, startperiod), frequency=freq );
ts.plot(ts3, gpars=list(ylab="Log[ LPI ]", col=c(1,2,3),frame = TRUE, axes=FALSE),
        main="Log of Originals and Trends");
abline(inter_o, slope_o, col="blue", lty=2);
abline(inter_o_nsn, slope_o_nsn, col="blue", lty=2);

legend(2005.0, 4.61,
       c("Original","Not SN jobs","Linear (LSQ) Trend fits", "X11 trend"),
       col=c("black","red","blue","green"),
       text.col="black",
       lty=c(1,1,2,1),
       lwd=c(1,1,1,1),
       pch=c(" ", " ", " ", " "),
       merge=FALSE,
       bty='n');

axis(1, axTicks(1), labels=xlabels, las=1)
axis(2);

#---------------------------------------------------------------------
# compute ratio of seasonal factors from difference of residuals and plot:

SBratio <- exp(resid_o_nsn - resid_o);

#bitmap(file="SBratios.png", width=5, height=5, res=400)
png(file="SBratios.png", width = 680, height = 450, pointsize = 12, bg="white", res = 200);

startyear <- 2003;
startperiod <- 2;
freq <- 4;
ts4 <- ts( SBratio, start=c(startyear, startperiod), frequency=freq );
ts4;
ts.plot(ts4, gpars=list(ylab="S(NSN) / S(total)", col=c(1), frame = TRUE, axes=FALSE),
        main="Seasonal factor ratio: \"non-safety net jobs / total\"");

axis(1, axTicks(1), labels=xlabels, las=1)
axis(2);

#---------------------------------------------------------------------
# Compute mvts in trends in % for use in method below:

t_o <- exp(o - resid_o);
t_o_nsn <- exp(o_nsn - resid_o_nsn);
mvt_trend_o <- 100*diff(t_o,1)/t_o[1:(length(t_o)-1)];
mvt_trend_o_nsn <- 100*diff(t_o_nsn,1)/t_o_nsn[1:(length(t_o_nsn)-1)];

#---------------------------------------------------------------------
# Estimate uncertainty in each trend movement in % and also difference
# |mvt_trend_o - mvt_trend_o_nsn|:

err_slope_o <- par_o[2,2];
err_slope_o_nsn <- par_o_nsn[2,2];

Dt <- t_year[2] - t_year[1];
err_mvt_trend_o <- 100*Dt*exp(slope_o*Dt)*err_slope_o;
err_mvt_trend_o_nsn <- 100*Dt*exp(slope_o_nsn*Dt)*err_slope_o_nsn;

diff_mvt_trends <- mvt_trend_o[2] - mvt_trend_o_nsn[1];
err_diff_mvt_trends <- sqrt(((err_mvt_trend_o)**2) + ((err_mvt_trend_o_nsn)**2));
tdiff <- diff_mvt_trends / err_diff_mvt_trends;
pvalue <- 2*pt(tdiff, (length(mvt_trend_o)-1), lower.tail = FALSE);

write(paste("\nDiff_trends =",diff_mvt_trends,"+/-",err_diff_mvt_trends,
            "\nPvalue =",100*pvalue,"%\n"),file="");
err_slope_o; err_slope_o_nsn;
#---------------------------------------------------------------------
# Compute potential trend break for periods at current end of series:

t_year_nosna <- 1996.75;  # => Q4 1996
t_year_tb <- 2006.0 + 0.25*cumsum(rep(1,4));
TB <- exp( (slope_o - slope_o_nsn)*(t_year_tb - t_year_nosna) );

write(paste("Trend Break =",TB,"at timepoint:",t_year_tb), file="");

#---------------------------------------------------------------------
# Compute NEW trend break for periods at current end of series:

t_year_mod <- 2003.0 + 0.25*cumsum(rep(1,16));
modelT_o <- exp(slope_o*t_year_mod + inter_o);
modelT_NSN <- modelT_o / TB;

write(paste("\n"));

t_year_mod;
modelT_o;
modelT_NSN;

TBnewSep <- modelT_o[[14]] / (modelT_o[[13]] + (modelT_NSN[[14]] - modelT_NSN[[13]]) );
write(paste("\nNEW Trend Break =",TBnewSep,"at timepoint: Sep 2006\n"), file="");

#---------------------------------------------------------------------
# Compute seasonal factor (SB) ratios using global least squares method
# to estimate four <S_nsn/S> ratios across all years, i.e. one for each qtr. 

R_NSN <- ( 1 + (MvtNSN/100) ) / ( 1+ (MvtTot/100) );
R_NSN_trend <- ( 1 + (mvt_trend_o_nsn/100) ) / ( 1+ (mvt_trend_o/100) );
rnsn <- log(R_NSN/R_NSN_trend);

Ac1 <- c(5, -2, 0, -3);
Ac2 <- c(-2, 4, -2, 0);
Ac3 <- c(0, -2, 5, -3);
Ac4 <- c(-3, 0, -3, 6);
A <- cbind(Ac1, Ac2, Ac3, Ac4);

s1 <- rnsn[[3]] + rnsn[[7]] + rnsn[[11]] - rnsn[[4]] - rnsn[[8]];
s2 <- rnsn[[4]] + rnsn[[8]] - rnsn[[5]] - rnsn[[9]];
s3 <- rnsn[[5]] + rnsn[[9]] - rnsn[[2]] - rnsn[[6]] - rnsn[[10]];
s4 <- -rnsn[[3]] - rnsn[[7]] - rnsn[[11]] + rnsn[[2]] + rnsn[[6]] + rnsn[[10]];
B <- c(s1, s2, s3, s4);

# Since above system [A.X = B] is singular (i.e. det(A)=0),
# use SVD to find closest solution vector:

svd <- svd(A);
#Asanitycheck <- svd$u %*% diag(svd$d) %*% t(svd$v);
invD <- diag(1/svd$d); # take inverse of diagonal matrix from svd.
invD[4,4] = 0;         # set the singular element to zero.
X <- svd$v %*% invD %*% t(svd$u) %*% B;

exp(X);

# Check that residuals [R = A.X - B ~ 0] appropriately small:
# A %*% X - B;