# SURexamples.ssc			script file for examples used in SUR chapter
#
# author: Eric Zivot
# created: December 18, 2001
# revised: March 31, 2001
#
# Examples: 
#
# 1. Exchange rate system
# 2. CAPM system
#
# Notes:
#
# 1. For SUR and NLSUR to work correctly, the data in the data
# slot of a timeSeries object must be a data frame. If the data
# are a matrix then the name information gets lost.
####################################################################

#
# Example:  SUR estimation of exchange rate system
#

formula.list = list(USCNS.diff~USCN.FP.lag1,
USDMS.diff~USDM.FP.lag1,
USFRS.diff~USFR.FP.lag1,
USILS.diff~USIL.FP.lag1,
USJYS.diff~USJY.FP.lag1,
USUKS.diff~USUK.FP.lag1)

# SUR estimation w/o iteration. Note: start sample in 1978 to eliminate NA values
# for USJY data

args(SUR)
sur.fit = SUR(formula.list,data=surex1.ts,
start="Aug 1978",in.format="%m %Y")
class(sur.fit)
names(sur.fit)
sur.fit
summary(sur.fit)

# compute iterated estimator

sur.fit2 = SUR(formula.list,data=surex1.ts,
start="Aug 1978",in.format="%m %Y",iterate=T)
coef(sur.fit)
coef(sur.fit2)
# alternatively, use
cbind(coef(sur.fit),coef(sur.fit2))

# illustrate plot method

# make sure trellis color options are set for white background
plot(sur.fit, which=4)
# select option 5 to compute residual ACF plot

#
# check residual correlation 
#
sd.vals = sqrt(diag(sur.fit$Sigma))
cor.mat = sur.fit$Sigma/outer(sd.vals,sd.vals)
cor.mat

#
# compare to OLS fits
#
ols.uscn.fit = OLS(USCNS.diff~USCN.FP.lag1,data=surex1.ts)
summary(ols.uscn.fit)

#
# compute Wald statistic.
#
bigR = matrix(0,6,12)
bigR[1,2] = bigR[2,4] = bigR[3,6] = bigR[4,8] = bigR[5,10] = bigR[6,12] = 1
rr = rep(1,6)
bHat = as.vector(coef(sur.fit))
avar = bigR%*%vcov(sur.fit)%*%t(bigR)
Wald = t((bigR%*%bHat-rr))%*%solve(avar)%*%(bigR%*%bHat-rr)
Wald
1-pchisq(Wald,6)

#
# compute LR statistic
#

# estimate restricted model
formula.list = list((USCNS.diff-USCN.FP.lag1)~1,
(USDMS.diff-USDM.FP.lag1)~1,
(USFRS.diff-USFR.FP.lag1)~1,
(USILS.diff-USIL.FP.lag1)~1,
(USJYS.diff-USJY.FP.lag1)~1,
(USUKS.diff-USUK.FP.lag1)~1)
sur.fit2r = SUR(formula.list,data=surex1.ts,
start="Aug 1978",in.format="%m %Y",iterate=T)
sur.fit2r

# compute LR statistic
nobs = nrow(residuals(sur.fit2r))
LR = nobs*(determinant(sur.fit2r$Sigma,log=T)$modulus-
determinant(sur.fit2$Sigma,log=T)$modulus)
as.numeric(LR)
1-pchisq(LR,6)

#
# Example: Estimating and Testing the CAPM
#

# Do LR test of CAPM by estimating restricted and unrestriced SUR model
# should be able to do this by specifying no constant in the regression with -1
# in formula.

colIds(black.ts)
formula.list = list(BOISE~MARKET,
CITCRP~MARKET,
CONED~MARKET,
CONTIL~MARKET,
DATGEN~MARKET,
DEC~MARKET,
DELTA~MARKET,
GENMIL~MARKET,
GERBER~MARKET,
IBM~MARKET,
MOBIL~MARKET,
PANAM~MARKET,
PSNH~MARKET,
TANDY~MARKET,
TEXACO~MARKET,
WEYER~MARKET)

capm.fit = SUR(formula.list,data=black.ts)
coef(capm.fit)

# construct Wald statistic to test CAPM

# fit restricted model for CAPM estimation
restricted.formula.list = list(BOISE~MARKET-1,
CITCRP~MARKET-1,
CONED~MARKET-1,
CONTIL~MARKET-1,
DATGEN~MARKET-1,
DEC~MARKET-1,
DELTA~MARKET-1,
GENMIL~MARKET-1,
GERBER~MARKET-1,
IBM~MARKET-1,
MOBIL~MARKET-1,
PANAM~MARKET-1,
PSNH~MARKET-1,
TANDY~MARKET-1,
TEXACO~MARKET-1,
WEYER~MARKET-1)

capm.restricted.fit = SUR(restricted.formula.list,data=black.ts,
iter=T)

nobs = nrow(residuals(capm.restricted.fit))
LR = nobs*(determinant(capm.fit$Sigma,log=T)$modulus-
determinant(capm.restricted.fit$Sigma,log=T)$modulus)
LR

#
# Nonlinear SUR models
#
args(NLSUR)
formula.list = list(y1~b10+b11*x1^b,
y2~b20+b21*x1^b,
y3~b30+b31*x1^b)
start.vals = c(0,1,0,1,0,1,0.5)
names(start.vals) = c("b10","b11","b20","b21","b30","b31","b")

test.dat = mk.zero1[,1:4]
formula.list = list(M.3~C1+C2*M.1^C3)
nlsur.fit = NLSUR(formula.list,data=test.dat)

# estimate Black form of CAPM using berndt data
colIds(black.ts)
# create formulas for restricted model
formula.list = list(BOISE~(1-b1)*g + b1*MARKET,
CITCRP~(1-b2)*g + b2*MARKET,
CONED~(1-b3)*g + b3*MARKET,
CONTIL~(1-b4)*g + b4*MARKET,
DATGEN~(1-b5)*g + b5*MARKET,
DEC~(1-b6)*g + b6*MARKET,
DELTA~(1-b7)*g + b7*MARKET,
GENMIL~(1-b8)*g + b8*MARKET,
GERBER~(1-b9)*g + b9*MARKET,
IBM~(1-b10)*g + b10*MARKET,
MOBIL~(1-b11)*g + b11*MARKET,
PANAM~(1-b12)*g + b12*MARKET,
PSNH~(1-b13)*g + b13*MARKET,
TANDY~(1-b14)*g + b14*MARKET,
TEXACO~(1-b15)*g + b15*MARKET,
WEYER~(1-b16)*g + b16*MARKET)

start.vals = c(0,rep(1,16))
names(start.vals) = c("g",paste("b",1:16,sep=""))

nlsur.fit = NLSUR(formula.list,data=black.ts,
coef=start.vals,start="Jan 1983",in.format="%m %Y")
names(nlsur.fit)
nlsur.fit$message
nlsur.fit
 
# compute standard error for b
std.ers = sqrt(diag(vcov(nlsur.fit)))
cbind(coef(nlsur.fit),std.ers)

plot(nlsur.fit, ask=F)

# unrestricted model
colIds(black.ts)
formula.list = list(BOISE~a1+b1*MARKET,
CITCRP~a2+b2*MARKET,
CONED~a3+b3*MARKET,
CONTIL~a4+b4*MARKET,
DATGEN~a5+b5*MARKET,
DEC~a6+b6*MARKET,
DELTA~a7+b7*MARKET,
GENMIL~a8+b8*MARKET,
GERBER~a9+b9*MARKET,
IBM~a10+b10*MARKET,
MOBIL~a11+b11*MARKET,
PANAM~a12+b12*MARKET,
PSNH~a13+b13*MARKET,
TANDY~a14+b14*MARKET,
TEXACO~a15+b15*MARKET,
WEYER~a16+b16*MARKET)

start.vals = c(rep(0,16),rep(1,16))
names(start.vals) = 
c(paste("a",1:16,sep=""),paste("b",1:16,sep=""))

nlsur.fit2 = NLSUR(formula.list,data=black.ts,
coef=start.vals,start="Jan 1983",in.format="%m %Y")

# compute LR statistic
nobs = nrow(residuals(nlsur.fit2))
LR = nobs*(determinant(nlsur.fit$Sigma,log=T)$modulus-
determinant(nlsur.fit2$Sigma,log=T)$modulus)
as.numeric(LR)
1-pchisq(LR,16)

#
# Example: Estimating and testing Factor model a la Burmiester and Elton
# see JBES 1988 article.
#

#
# estimation of exchange rate system imposing cross equation restrictions
#
# NLSUR estimation of exchange rate system subject to cross equation restrictions
# 
formula.list = list(USCNS.diff~a1+g*USCN.FP.lag1,
USDMS.diff~a2+g*USDM.FP.lag1,
USFRS.diff~a3+g*USFR.FP.lag1,
USILS.diff~a4+g*USIL.FP.lag1,
USJYS.diff~a5+g*USJY.FP.lag1,
USUKS.diff~a6+g*USUK.FP.lag1)

start.vals = c(rep(0,6),1)
names(start.vals) = c(paste("a",1:6,sep=""),"g")

# NLSUR estimation. Note: start sample in 1978 to eliminate NA values
# for USJY data
args(NLSUR)
nlsur.fit = NLSUR(formula.list,data=surex1.ts,
start="Aug 1978",in.format="%m %Y")
class(nlsur.fit)
names(nlsur.fit)
nlsur.fit

# NLSUR estimation specifying starting values
nlsur.fit = NLSUR(formula.list,data=surex1.ts,coef=start.vals,
start="Aug 1978",in.format="%m %Y")
nlsur.fit

# compute standard error for g
sqrt(diag(vcov(nlsur.fit)))[7]

# LR test of common parameter restriction
nobs = nrow(residuals(nlsur.fit))
LR = nobs*(determinant(nlsur.fit$Sigma,log=T)$modulus
-determinant(sur.fit2$Sigma,log=T)$modulus)
as.numeric(LR)
1 - pchisq(LR,6)