x <- array(list(235.1
,280.7
,264.6
,240.7
,201.4
,240.8
,241.1
,223.8
,206.1
,174.7
,203.3
,220.5
,299.5
,347.4
,338.3
,327.7
,351.6
,396.6
,438.8
,395.6
,363.5
,378.8
,357
,369
,464.8
,479.1
,431.3
,366.5
,326.3
,355.1
,331.6
,261.3
,249
,205.5
,235.6
,240.9
,264.9
,253.8
,232.3
,193.8
,177
,213.2
,207.2
,180.6
,188.6
,175.4
,199
,179.6
,225.8
,234
,200.2
,183.6
,178.2
,203.2
,208.5
,191.8
,172.8
,148
,159.4
,154.5
,213.2
,196.4
,182.8
,176.4
,153.6
,173.2
,171
,151.2
,161.9
,157.2
,201.7
,236.4
,356.1
,398.3
,403.7
,384.6
,365.8
,368.1
,367.9
,347
,343.3
,292.9
,311.5
,300.9
,366.9
,356.9
,329.7
,316.2
,269
,289.3
,266.2
,253.6
,233.8
,228.4
,253.6
,260.1
,306.6
,309.2
,309.5
,271
,279.9
,317.9
,298.4
,246.7
,227.3
,209.1
,259.9
,266
,320.6
,308.5
,282.2
,262.7
,263.5
,313.1
,284.3
,252.6
,250.3
,246.5
,312.7
,333.2
,446.4
,511.6
,515.5
,506.4
,483.2
,522.3
,509.8
,460.7
,405.8
,375
,378.5
,406.8
,467.8
,469.8
,429.8
,355.8
,332.7
,378
,360.5
,334.7
,319.5
,323.1
,363.6
,352.1
,411.9
,388.6
,416.4
,360.7
,338
,417.2
,388.4
,371.1
,331.5
,353.7
,396.7
,447
,533.5
,565.4
,542.3
,488.7
,467.1
,531.3
,496.1
,444
,403.4
,386.3
,394.1
,404.1
,462.1
,448.1
,432.3
,386.3
,395.2
,421.9
,382.9
,384.2
,345.5
,323.4
,372.6
,376
,462.7
,487
,444.2
,399.3
,394.9
,455.4
,414
,375.5
,347
,339.4
,385.8
,378.8
,451.8
,446.1
,422.5
,383.1
,352.8
,445.3
,367.5
,355.1
,326.2
,319.8
,331.8
,340.9
,394.1
,417.2
,369.9
,349.2
,321.4
,405.7
,342.9
,316.5
,284.2
,270.9
,288.8
,278.8
,324.4
,310.9
,299
,273
,279.3
,359.2
,305
,282.1
,250.3
,246.5
,257.9
,266.5
,315.9
,318.4
,295.4
,266.4
,245.8
,362.8
,324.9
,294.2
,289.5
,295.2
,290.3
,272
,307.4
,328.7
,292.9
,249.1
,230.4
,361.5
,321.7
,277.2
,260.7
,251
,257.6
,241.8
,287.5
,292.3
,274.7
,254.2
,230
,339
,318.2
,287
,295.8
,284
,271
,262.7
,340.6
,379.4
,373.3
,355.2
,338.4
,466.9
,451
,422
,429.2
,425.9
,460.7
,463.6
,541.4
,544.2
,517.5
,469.4
,439.4
,549
,533
,506.1
,484
,457
,481.5
,469.5
,544.7
,541.2
,521.5
,469.7
,434.4
,542.6
,517.3
,485.7
,465.8
,447
,426.6
,411.6
,467.5
,484.5
,451.2
,417.4
,379.9
,484.7
,455
,420.8
,416.5
,376.3
,405.6
,405.8
,500.8
,514
,475.5
,430.1
,414.4
,538
,526
,488.5
,520.2
,504.4
,568.5
,610.6
,818
,830.9
,835.9
,782
,762.3
,856.9
,820.9
,769.6
,752.2
,724.4
,723.1
,719.5
,817.4
,803.3
,752.5
,689
,630.4
,765.5
,757.7
,732.2
,702.6
,683.3
,709.5
,702.2
,784.8
,810.9
,755.6
,656.8
,615.1
,745.3
,694.1
,675.7
,643.7
,622.1
,634.6
,588
,689.7
,673.9
,647.9
,568.8
,545.7
,632.6
,643.8
,593.1
,579.7
,546
,562.9
,572.5)
,dim=c(1
,372)
,dimnames=list(c('Unemployment')
,1:372))
 y <- array(NA,dim=c(1,372),dimnames=list(c('Unemployment'),1:372))
 for (i in 1:dim(x)[1])
 {
 	for (j in 1:dim(x)[2])
 	{
 		y[i,j] <- as.numeric(x[i,j])
 	}
 }
par3 = 'Linear Trend'
par2 = 'Include Monthly Dummies'
par1 = '1'
library(lattice)
library(lmtest)
n25 <- 25 #minimum number of obs. for Goldfeld-Quandt test
par1 <- as.numeric(par1)
x <- t(y)
k <- length(x[1,])
n <- length(x[,1])
x1 <- cbind(x[,par1], x[,1:k!=par1])
mycolnames <- c(colnames(x)[par1], colnames(x)[1:k!=par1])
colnames(x1) <- mycolnames #colnames(x)[par1]
x <- x1
if (par3 == 'First Differences'){
x2 <- array(0, dim=c(n-1,k), dimnames=list(1:(n-1), paste('(1-B)',colnames(x),sep='')))
for (i in 1:n-1) {
for (j in 1:k) {
x2[i,j] <- x[i+1,j] - x[i,j]
}
}
x <- x2
}
if (par2 == 'Include Monthly Dummies'){
x2 <- array(0, dim=c(n,11), dimnames=list(1:n, paste('M', seq(1:11), sep ='')))
for (i in 1:11){
x2[seq(i,n,12),i] <- 1
}
x <- cbind(x, x2)
}
if (par2 == 'Include Quarterly Dummies'){
x2 <- array(0, dim=c(n,3), dimnames=list(1:n, paste('Q', seq(1:3), sep ='')))
for (i in 1:3){
x2[seq(i,n,4),i] <- 1
}
x <- cbind(x, x2)
}
k <- length(x[1,])
if (par3 == 'Linear Trend'){
x <- cbind(x, c(1:n))
colnames(x)[k+1] <- 't'
}
x
k <- length(x[1,])
df <- as.data.frame(x)
(mylm <- lm(df))
(mysum <- summary(mylm))
if (n > n25) {
kp3 <- k + 3
nmkm3 <- n - k - 3
gqarr <- array(NA, dim=c(nmkm3-kp3+1,3))
numgqtests <- 0
numsignificant1 <- 0
numsignificant5 <- 0
numsignificant10 <- 0
for (mypoint in kp3:nmkm3) {
j <- 0
numgqtests <- numgqtests + 1
for (myalt in c('greater', 'two.sided', 'less')) {
j <- j + 1
gqarr[mypoint-kp3+1,j] <- gqtest(mylm, point=mypoint, alternative=myalt)$p.value
}
if (gqarr[mypoint-kp3+1,2] < 0.01) numsignificant1 <- numsignificant1 + 1
if (gqarr[mypoint-kp3+1,2] < 0.05) numsignificant5 <- numsignificant5 + 1
if (gqarr[mypoint-kp3+1,2] < 0.10) numsignificant10 <- numsignificant10 + 1
}
gqarr
}
postscript(file="/var/wessaorg/rcomp/tmp/1il3c1322584887.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) 
plot(x[,1], type='l', main='Actuals and Interpolation', ylab='value of Actuals and Interpolation (dots)', xlab='time or index')
points(x[,1]-mysum$resid)
grid()
dev.off()
postscript(file="/var/wessaorg/rcomp/tmp/2zy0r1322584887.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) 
plot(mysum$resid, type='b', pch=19, main='Residuals', ylab='value of Residuals', xlab='time or index')
grid()
dev.off()
postscript(file="/var/wessaorg/rcomp/tmp/3h93s1322584887.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) 
hist(mysum$resid, main='Residual Histogram', xlab='values of Residuals')
grid()
dev.off()
postscript(file="/var/wessaorg/rcomp/tmp/44q2s1322584887.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) 
densityplot(~mysum$resid,col='black',main='Residual Density Plot', xlab='values of Residuals')
dev.off()
postscript(file="/var/wessaorg/rcomp/tmp/5l3e61322584887.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) 
qqnorm(mysum$resid, main='Residual Normal Q-Q Plot')
qqline(mysum$resid)
grid()
dev.off()
(myerror <- as.ts(mysum$resid))
postscript(file="/var/wessaorg/rcomp/tmp/6omvj1322584887.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) 
dum <- cbind(lag(myerror,k=1),myerror)
dum
dum1 <- dum[2:length(myerror),]
dum1
z <- as.data.frame(dum1)
z
plot(z,main=paste('Residual Lag plot, lowess, and regression line'), ylab='values of Residuals', xlab='lagged values of Residuals')
lines(lowess(z))
abline(lm(z))
grid()
dev.off()
postscript(file="/var/wessaorg/rcomp/tmp/7bbsd1322584887.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) 
acf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Autocorrelation Function')
grid()
dev.off()
postscript(file="/var/wessaorg/rcomp/tmp/884zr1322584887.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) 
pacf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Partial Autocorrelation Function')
grid()
dev.off()
postscript(file="/var/wessaorg/rcomp/tmp/9q34q1322584887.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) 
opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
plot(mylm, las = 1, sub='Residual Diagnostics')
par(opar)
dev.off()
if (n > n25) {
postscript(file="/var/wessaorg/rcomp/tmp/1014ij1322584887.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) 
plot(kp3:nmkm3,gqarr[,2], main='Goldfeld-Quandt test',ylab='2-sided p-value',xlab='breakpoint')
grid()
dev.off()
}

#Note: the /var/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
load(file="/var/wessaorg/rcomp/createtable")

a<-table.start()
a<-table.row.start(a)
a<-table.element(a, 'Multiple Linear Regression - Estimated Regression Equation', 1, TRUE)
a<-table.row.end(a)
myeq <- colnames(x)[1]
myeq <- paste(myeq, '[t] = ', sep='')
for (i in 1:k){
if (mysum$coefficients[i,1] > 0) myeq <- paste(myeq, '+', '')
myeq <- paste(myeq, mysum$coefficients[i,1], sep=' ')
if (rownames(mysum$coefficients)[i] != '(Intercept)') {
myeq <- paste(myeq, rownames(mysum$coefficients)[i], sep='')
if (rownames(mysum$coefficients)[i] != 't') myeq <- paste(myeq, '[t]', sep='')
}
}
myeq <- paste(myeq, ' + e[t]')
a<-table.row.start(a)
a<-table.element(a, myeq)
a<-table.row.end(a)
a<-table.end(a)
table.save(a,file="/var/wessaorg/rcomp/tmp/11dlci1322584887.tab") 
a<-table.start()
a<-table.row.start(a)
a<-table.element(a,hyperlink('http://www.xycoon.com/ols1.htm','Multiple Linear Regression - Ordinary Least Squares',''), 6, TRUE)
a<-table.row.end(a)
a<-table.row.start(a)
a<-table.element(a,'Variable',header=TRUE)
a<-table.element(a,'Parameter',header=TRUE)
a<-table.element(a,'S.D.',header=TRUE)
a<-table.element(a,'T-STAT<br />H0: parameter = 0',header=TRUE)
a<-table.element(a,'2-tail p-value',header=TRUE)
a<-table.element(a,'1-tail p-value',header=TRUE)
a<-table.row.end(a)
for (i in 1:k){
a<-table.row.start(a)
a<-table.element(a,rownames(mysum$coefficients)[i],header=TRUE)
a<-table.element(a,mysum$coefficients[i,1])
a<-table.element(a, round(mysum$coefficients[i,2],6))
a<-table.element(a, round(mysum$coefficients[i,3],4))
a<-table.element(a, round(mysum$coefficients[i,4],6))
a<-table.element(a, round(mysum$coefficients[i,4]/2,6))
a<-table.row.end(a)
}
a<-table.end(a)
table.save(a,file="/var/wessaorg/rcomp/tmp/12ceax1322584887.tab") 
a<-table.start()
a<-table.row.start(a)
a<-table.element(a, 'Multiple Linear Regression - Regression Statistics', 2, TRUE)
a<-table.row.end(a)
a<-table.row.start(a)
a<-table.element(a, 'Multiple R',1,TRUE)
a<-table.element(a, sqrt(mysum$r.squared))
a<-table.row.end(a)
a<-table.row.start(a)
a<-table.element(a, 'R-squared',1,TRUE)
a<-table.element(a, mysum$r.squared)
a<-table.row.end(a)
a<-table.row.start(a)
a<-table.element(a, 'Adjusted R-squared',1,TRUE)
a<-table.element(a, mysum$adj.r.squared)
a<-table.row.end(a)
a<-table.row.start(a)
a<-table.element(a, 'F-TEST (value)',1,TRUE)
a<-table.element(a, mysum$fstatistic[1])
a<-table.row.end(a)
a<-table.row.start(a)
a<-table.element(a, 'F-TEST (DF numerator)',1,TRUE)
a<-table.element(a, mysum$fstatistic[2])
a<-table.row.end(a)
a<-table.row.start(a)
a<-table.element(a, 'F-TEST (DF denominator)',1,TRUE)
a<-table.element(a, mysum$fstatistic[3])
a<-table.row.end(a)
a<-table.row.start(a)
a<-table.element(a, 'p-value',1,TRUE)
a<-table.element(a, 1-pf(mysum$fstatistic[1],mysum$fstatistic[2],mysum$fstatistic[3]))
a<-table.row.end(a)
a<-table.row.start(a)
a<-table.element(a, 'Multiple Linear Regression - Residual Statistics', 2, TRUE)
a<-table.row.end(a)
a<-table.row.start(a)
a<-table.element(a, 'Residual Standard Deviation',1,TRUE)
a<-table.element(a, mysum$sigma)
a<-table.row.end(a)
a<-table.row.start(a)
a<-table.element(a, 'Sum Squared Residuals',1,TRUE)
a<-table.element(a, sum(myerror*myerror))
a<-table.row.end(a)
a<-table.end(a)
table.save(a,file="/var/wessaorg/rcomp/tmp/13hp461322584887.tab") 
a<-table.start()
a<-table.row.start(a)
a<-table.element(a, 'Multiple Linear Regression - Actuals, Interpolation, and Residuals', 4, TRUE)
a<-table.row.end(a)
a<-table.row.start(a)
a<-table.element(a, 'Time or Index', 1, TRUE)
a<-table.element(a, 'Actuals', 1, TRUE)
a<-table.element(a, 'Interpolation<br />Forecast', 1, TRUE)
a<-table.element(a, 'Residuals<br />Prediction Error', 1, TRUE)
a<-table.row.end(a)
for (i in 1:n) {
a<-table.row.start(a)
a<-table.element(a,i, 1, TRUE)
a<-table.element(a,x[i])
a<-table.element(a,x[i]-mysum$resid[i])
a<-table.element(a,mysum$resid[i])
a<-table.row.end(a)
}
a<-table.end(a)
table.save(a,file="/var/wessaorg/rcomp/tmp/14czpn1322584887.tab") 
if (n > n25) {
a<-table.start()
a<-table.row.start(a)
a<-table.element(a,'Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
a<-table.row.end(a)
a<-table.row.start(a)
a<-table.element(a,'p-values',header=TRUE)
a<-table.element(a,'Alternative Hypothesis',3,header=TRUE)
a<-table.row.end(a)
a<-table.row.start(a)
a<-table.element(a,'breakpoint index',header=TRUE)
a<-table.element(a,'greater',header=TRUE)
a<-table.element(a,'2-sided',header=TRUE)
a<-table.element(a,'less',header=TRUE)
a<-table.row.end(a)
for (mypoint in kp3:nmkm3) {
a<-table.row.start(a)
a<-table.element(a,mypoint,header=TRUE)
a<-table.element(a,gqarr[mypoint-kp3+1,1])
a<-table.element(a,gqarr[mypoint-kp3+1,2])
a<-table.element(a,gqarr[mypoint-kp3+1,3])
a<-table.row.end(a)
}
a<-table.end(a)
table.save(a,file="/var/wessaorg/rcomp/tmp/15bi0k1322584887.tab") 
a<-table.start()
a<-table.row.start(a)
a<-table.element(a,'Meta Analysis of Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
a<-table.row.end(a)
a<-table.row.start(a)
a<-table.element(a,'Description',header=TRUE)
a<-table.element(a,'# significant tests',header=TRUE)
a<-table.element(a,'% significant tests',header=TRUE)
a<-table.element(a,'OK/NOK',header=TRUE)
a<-table.row.end(a)
a<-table.row.start(a)
a<-table.element(a,'1% type I error level',header=TRUE)
a<-table.element(a,numsignificant1)
a<-table.element(a,numsignificant1/numgqtests)
if (numsignificant1/numgqtests < 0.01) dum <- 'OK' else dum <- 'NOK'
a<-table.element(a,dum)
a<-table.row.end(a)
a<-table.row.start(a)
a<-table.element(a,'5% type I error level',header=TRUE)
a<-table.element(a,numsignificant5)
a<-table.element(a,numsignificant5/numgqtests)
if (numsignificant5/numgqtests < 0.05) dum <- 'OK' else dum <- 'NOK'
a<-table.element(a,dum)
a<-table.row.end(a)
a<-table.row.start(a)
a<-table.element(a,'10% type I error level',header=TRUE)
a<-table.element(a,numsignificant10)
a<-table.element(a,numsignificant10/numgqtests)
if (numsignificant10/numgqtests < 0.1) dum <- 'OK' else dum <- 'NOK'
a<-table.element(a,dum)
a<-table.row.end(a)
a<-table.end(a)
table.save(a,file="/var/wessaorg/rcomp/tmp/16n0cd1322584887.tab") 
}

try(system("convert tmp/1il3c1322584887.ps tmp/1il3c1322584887.png",intern=TRUE))
try(system("convert tmp/2zy0r1322584887.ps tmp/2zy0r1322584887.png",intern=TRUE))
try(system("convert tmp/3h93s1322584887.ps tmp/3h93s1322584887.png",intern=TRUE))
try(system("convert tmp/44q2s1322584887.ps tmp/44q2s1322584887.png",intern=TRUE))
try(system("convert tmp/5l3e61322584887.ps tmp/5l3e61322584887.png",intern=TRUE))
try(system("convert tmp/6omvj1322584887.ps tmp/6omvj1322584887.png",intern=TRUE))
try(system("convert tmp/7bbsd1322584887.ps tmp/7bbsd1322584887.png",intern=TRUE))
try(system("convert tmp/884zr1322584887.ps tmp/884zr1322584887.png",intern=TRUE))
try(system("convert tmp/9q34q1322584887.ps tmp/9q34q1322584887.png",intern=TRUE))
try(system("convert tmp/1014ij1322584887.ps tmp/1014ij1322584887.png",intern=TRUE))