library(MASS)

ps.options(horizontal=FALSE,onefile=FALSE)
ps.options(pagecentre=TRUE)
ps.options(paper="letter")
ps.options(width=7.0,height=7.0)

snpname  = "snp.big_down.plt"
npbname  = "npb.big_down.plt"
plotname = "crash_rally.eps"

M <- 2
ngrd <- 50
l <- (2*ngrd+1)
lplt <- (M+1)*l^M

postscript(file=plotname)
par(mfrow=c(2,2),mar=c(2.5,4,1.5,2)+0.1) # mar=c(b,l,t,r)

plt  <- scan(snpname)
mu   <- plt[seq(from=1,to=M)]
sig  <- plt[seq(from=M+1,to=M+M*M)]
zinc <- plt[seq(from=M+M*M+1,to=M+M*M+M)]
plt  <- plt[seq(from=M+M*M+M+1,to=M+M*M+M+lplt)]

h <- matrix(plt[seq(from=M+1,to=lplt,by=M+1)],nrow=l,ncol=l)
snpx <- plt[seq(from=1,to=l*(M+1),by=M+1)]     
snpy <- plt[seq(from=2,to=lplt,by=l*(M+1))]

one <- seq(1,1,length=l)

snphx <- zinc[2]*h%*%one
snphy <- zinc[1]*one%*%h

plt  <- scan(npbname)
mu   <- plt[seq(from=1,to=M)]
sig  <- plt[seq(from=M+1,to=M+M*M)]
zinc <- plt[seq(from=M+M*M+1,to=M+M*M+M)]
plt  <- plt[seq(from=M+M*M+M+1,to=M+M*M+M+lplt)]

h <- matrix(plt[seq(from=M+1,to=lplt,by=M+1)],nrow=l,ncol=l)
npbx <- plt[seq(from=1,to=l*(M+1),by=M+1)]     
npby <- plt[seq(from=2,to=lplt,by=l*(M+1))]

one <- seq(1,1,length=l)

npbhx <- zinc[2]*h%*%one
npbhy <- zinc[1]*one%*%h

left <- min(c(snpx,npbx,snpy,npby))
rite <- max(c(snpx,npbx,snpy,npby))
top  <- max(c(snphx,npbhx,snphy,npbhy))

plot(c(left,rite),c(0,top),type='n',xlab="",ylab="") # y_1
lines(snpx,snphx,lty="dashed")
lines(npbx,npbhx,lty="solid")
title("S&P500 after crash")

plot(c(left,rite),c(0,top),type='n',xlab="",ylab="") # y_2
lines(snpy,snphy,lty="dashed")
lines(npby,npbhy,lty="solid")
title("NASDAQ after crash")

snpname  = "snp.big_up.plt"
npbname  = "npb.big_up.plt"

M <- 2
ngrd <- 50
l <- (2*ngrd+1)
lplt <- (M+1)*l^M

plt  <- scan(snpname)
mu   <- plt[seq(from=1,to=M)]
sig  <- plt[seq(from=M+1,to=M+M*M)]
zinc <- plt[seq(from=M+M*M+1,to=M+M*M+M)]
plt  <- plt[seq(from=M+M*M+M+1,to=M+M*M+M+lplt)]

h <- matrix(plt[seq(from=M+1,to=lplt,by=M+1)],nrow=l,ncol=l)
snpx <- plt[seq(from=1,to=l*(M+1),by=M+1)]     
snpy <- plt[seq(from=2,to=lplt,by=l*(M+1))]

one <- seq(1,1,length=l)

snphx <- zinc[2]*h%*%one
snphy <- zinc[1]*one%*%h

plt  <- scan(npbname)
mu   <- plt[seq(from=1,to=M)]
sig  <- plt[seq(from=M+1,to=M+M*M)]
zinc <- plt[seq(from=M+M*M+1,to=M+M*M+M)]
plt  <- plt[seq(from=M+M*M+M+1,to=M+M*M+M+lplt)]

h <- matrix(plt[seq(from=M+1,to=lplt,by=M+1)],nrow=l,ncol=l)
npbx <- plt[seq(from=1,to=l*(M+1),by=M+1)]     
npby <- plt[seq(from=2,to=lplt,by=l*(M+1))]

one <- seq(1,1,length=l)

npbhx <- zinc[2]*h%*%one
npbhy <- zinc[1]*one%*%h

left <- min(c(snpx,npbx,snpy,npby))
rite <- max(c(snpx,npbx,snpy,npby))
top  <- max(c(snphx,npbhx,snphy,npbhy))

plot(c(left,rite),c(0,top),type='n',xlab="",ylab="") # y_1
lines(snpx,snphx,lty="dashed")
lines(npbx,npbhx,lty="solid")
title("S&P500 after rally")

plot(c(left,rite),c(0,top),type='n',xlab="",ylab="") # y_2
lines(snpy,snphy,lty="dashed")
lines(npby,npbhy,lty="solid")
title("NASDAQ after rally")

dev.off()

