R version 2.12.0 (2010-10-15)
Copyright (C) 2010 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i486-pc-linux-gnu (32-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> x <- array(list(-999.0
+ ,38.6
+ ,6654.000
+ ,5712.000
+ ,645.0
+ ,3
+ ,5
+ ,3
+ ,6.3
+ ,4.5
+ ,1.000
+ ,6.600
+ ,42.0
+ ,3
+ ,1
+ ,3
+ ,-999.0
+ ,14.0
+ ,3.385
+ ,44.500
+ ,60.0
+ ,1
+ ,1
+ ,1
+ ,-999.0
+ ,-999.0
+ ,0.920
+ ,5.700
+ ,25.0
+ ,5
+ ,2
+ ,3
+ ,2.1
+ ,69.0
+ ,2547.000
+ ,4603.000
+ ,624.0
+ ,3
+ ,5
+ ,4
+ ,9.1
+ ,27.0
+ ,10.550
+ ,179.500
+ ,180.0
+ ,4
+ ,4
+ ,4
+ ,15.8
+ ,19.0
+ ,0.023
+ ,0.300
+ ,35.0
+ ,1
+ ,1
+ ,1
+ ,5.2
+ ,30.4
+ ,160.000
+ ,169.000
+ ,392.0
+ ,4
+ ,5
+ ,4
+ ,10.9
+ ,28.0
+ ,3.300
+ ,25.600
+ ,63.0
+ ,1
+ ,2
+ ,1
+ ,8.3
+ ,50.0
+ ,52.160
+ ,440.000
+ ,230.0
+ ,1
+ ,1
+ ,1
+ ,11.0
+ ,7.0
+ ,0.425
+ ,6.400
+ ,112.0
+ ,5
+ ,4
+ ,4
+ ,3.2
+ ,30.0
+ ,465.000
+ ,423.000
+ ,281.0
+ ,5
+ ,5
+ ,5
+ ,7.6
+ ,-999.0
+ ,0.550
+ ,2.400
+ ,-999.0
+ ,2
+ ,1
+ ,2
+ ,-999.0
+ ,40.0
+ ,187.100
+ ,419.000
+ ,365.0
+ ,5
+ ,5
+ ,5
+ ,6.3
+ ,0.075
+ ,1.200
+ ,42.0
+ ,1
+ ,1
+ ,1
+ ,8.6
+ ,50.0
+ ,3.000
+ ,25.000
+ ,28.0
+ ,2
+ ,2
+ ,2
+ ,6.6
+ ,6.0
+ ,0.785
+ ,3.500
+ ,42.0
+ ,2
+ ,2
+ ,2
+ ,9.5
+ ,10.4
+ ,0.200
+ ,5.000
+ ,120.0
+ ,2
+ ,2
+ ,2
+ ,4.8
+ ,34.0
+ ,1.410
+ ,17.500
+ ,-999.0
+ ,1
+ ,2
+ ,1
+ ,12.0
+ ,7.0
+ ,60.000
+ ,81.000
+ ,-999.0
+ ,1
+ ,1
+ ,1
+ ,-999.0
+ ,28.0
+ ,529.000
+ ,680.000
+ ,400.0
+ ,5
+ ,5
+ ,5
+ ,3.3
+ ,20.0
+ ,27.660
+ ,115.000
+ ,148.0
+ ,5
+ ,5
+ ,5
+ ,11.0
+ ,3.9
+ ,0.120
+ ,1.000
+ ,16.0
+ ,3
+ ,1
+ ,2
+ ,-999.0
+ ,39.3
+ ,207.000
+ ,406.000
+ ,252.0
+ ,1
+ ,4
+ ,1
+ ,4.7
+ ,41.0
+ ,85.000
+ ,325.000
+ ,310.0
+ ,1
+ ,3
+ ,1
+ ,-999.0
+ ,16.2
+ ,36.330
+ ,119.500
+ ,63.0
+ ,1
+ ,1
+ ,1
+ ,10.4
+ ,9.0
+ ,0.101
+ ,4.000
+ ,28.0
+ ,5
+ ,1
+ ,3
+ ,7.4
+ ,7.6
+ ,1.040
+ ,5.500
+ ,68.0
+ ,5
+ ,3
+ ,4
+ ,2.1
+ ,46.0
+ ,521.000
+ ,655.000
+ ,336.0
+ ,5
+ ,5
+ ,5
+ ,-999.0
+ ,22.4
+ ,100.000
+ ,157.000
+ ,100.0
+ ,1
+ ,1
+ ,1
+ ,-999.0
+ ,16.3
+ ,35.000
+ ,56.000
+ ,33.0
+ ,3
+ ,5
+ ,4
+ ,7.7
+ ,2.6
+ ,0.005
+ ,0.140
+ ,21.5
+ ,5
+ ,2
+ ,4
+ ,17.9
+ ,24.0
+ ,0.010
+ ,0.250
+ ,50.0
+ ,1
+ ,1
+ ,1
+ ,6.1
+ ,100.0
+ ,62.000
+ ,1320.000
+ ,267.0
+ ,1
+ ,1
+ ,1
+ ,8.2
+ ,-999.0
+ ,0.122
+ ,3.000
+ ,30.0
+ ,2
+ ,1
+ ,1
+ ,8.4
+ ,-999.0
+ ,1.350
+ ,8.100
+ ,45.0
+ ,3
+ ,1
+ ,3
+ ,11.9
+ ,3.2
+ ,0.023
+ ,0.400
+ ,19.0
+ ,4
+ ,1
+ ,3
+ ,10.8
+ ,2.0
+ ,0.048
+ ,0.330
+ ,30.0
+ ,4
+ ,1
+ ,3
+ ,13.8
+ ,5.0
+ ,1.700
+ ,6.300
+ ,12.0
+ ,2
+ ,1
+ ,1
+ ,14.3
+ ,6.5
+ ,3.500
+ ,10.800
+ ,120.0
+ ,2
+ ,1
+ ,1
+ ,-999.0
+ ,23.6
+ ,250.000
+ ,490.000
+ ,440.0
+ ,5
+ ,5
+ ,5
+ ,15.2
+ ,12.0
+ ,0.480
+ ,15.500
+ ,140.0
+ ,2
+ ,2
+ ,2
+ ,10.0
+ ,20.2
+ ,10.000
+ ,115.000
+ ,170.0
+ ,4
+ ,4
+ ,4
+ ,11.9
+ ,13.0
+ ,1.620
+ ,11.400
+ ,17.0
+ ,2
+ ,1
+ ,2
+ ,6.5
+ ,27.0
+ ,192.000
+ ,180.000
+ ,115.0
+ ,4
+ ,4
+ ,4
+ ,7.5
+ ,18.0
+ ,2.500
+ ,12.100
+ ,31.0
+ ,5
+ ,5
+ ,5
+ ,-999.0
+ ,13.7
+ ,4.288
+ ,39.200
+ ,63.0
+ ,2
+ ,2
+ ,2
+ ,10.6
+ ,4.7
+ ,0.280
+ ,1.900
+ ,21.0
+ ,3
+ ,1
+ ,3
+ ,7.4
+ ,9.8
+ ,4.235
+ ,50.400
+ ,52.0
+ ,1
+ ,1
+ ,1
+ ,8.4
+ ,29.0
+ ,6.800
+ ,179.000
+ ,164.0
+ ,2
+ ,3
+ ,2
+ ,5.7
+ ,7.0
+ ,0.750
+ ,12.300
+ ,225.0
+ ,2
+ ,2
+ ,2
+ ,4.9
+ ,6.0
+ ,3.600
+ ,21.000
+ ,225.0
+ ,3
+ ,2
+ ,3
+ ,-999.0
+ ,17.0
+ ,14.830
+ ,98.200
+ ,150.0
+ ,5
+ ,5
+ ,5
+ ,3.2
+ ,20.0
+ ,55.500
+ ,175.000
+ ,151.0
+ ,5
+ ,5
+ ,5
+ ,-999.0
+ ,12.7
+ ,1.400
+ ,12.500
+ ,90.0
+ ,2
+ ,2
+ ,2
+ ,8.1
+ ,3.5
+ ,0.060
+ ,1.000
+ ,-999.0
+ ,3
+ ,1
+ ,1
+ ,11.0
+ ,4.5
+ ,0.900
+ ,2.600
+ ,38.0
+ ,2
+ ,1
+ ,2
+ ,-999.0
+ ,7.5
+ ,2.000
+ ,17.000
+ ,200.0
+ ,3
+ ,1
+ ,3
+ ,13.2
+ ,2.3
+ ,0.104
+ ,2.500
+ ,46.0
+ ,3
+ ,2
+ ,2
+ ,9.7
+ ,13.0
+ ,4.050
+ ,58.000
+ ,210.0
+ ,4
+ ,3
+ ,4
+ ,12.8
+ ,3.0
+ ,3.500
+ ,3.900
+ ,14.0
+ ,2
+ ,1
+ ,1
+ ,4.9
+ ,24.0
+ ,4.190
+ ,12.300
+ ,60.0
+ ,3
+ ,1
+ ,2)
+ ,dim=c(8
+ ,62)
+ ,dimnames=list(c('SWS'
+ ,'L'
+ ,'wb'
+ ,'wbr'
+ ,'tg'
+ ,'P'
+ ,'S'
+ ,'D
')
+ ,1:62))
> y <- array(NA,dim=c(8,62),dimnames=list(c('SWS','L','wb','wbr','tg','P','S','D
'),1:62))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'No Linear Trend'
> par2 = 'Do not include Seasonal Dummies'
> par1 = '1'
> #'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
> #Author: Prof. Dr. P. Wessa
> #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #Technical description: Write here your technical program description (don't use hard returns!)
> library(lattice)
> library(lmtest)
Loading required package: zoo
> 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
SWS L wb wbr tg P S D\r
1 -999.0 38.600 6654.000 5712.0 645 3 5 3.0
2 6.3 4.500 1.000 6.6 42 3 1 3.0
3 -999.0 14.000 3.385 44.5 60 1 1 1.0
4 -999.0 -999.000 0.920 5.7 25 5 2 3.0
5 2.1 69.000 2547.000 4603.0 624 3 5 4.0
6 9.1 27.000 10.550 179.5 180 4 4 4.0
7 15.8 19.000 0.023 0.3 35 1 1 1.0
8 5.2 30.400 160.000 169.0 392 4 5 4.0
9 10.9 28.000 3.300 25.6 63 1 2 1.0
10 8.3 50.000 52.160 440.0 230 1 1 1.0
11 11.0 7.000 0.425 6.4 112 5 4 4.0
12 3.2 30.000 465.000 423.0 281 5 5 5.0
13 7.6 -999.000 0.550 2.4 -999 2 1 2.0
14 -999.0 40.000 187.100 419.0 365 5 5 5.0
15 6.3 0.075 1.200 42.0 1 1 1 8.6
16 50.0 3.000 25.000 28.0 2 2 2 6.6
17 6.0 0.785 3.500 42.0 2 2 2 9.5
18 10.4 0.200 5.000 120.0 2 2 2 4.8
19 34.0 1.410 17.500 -999.0 1 2 1 12.0
20 7.0 60.000 81.000 -999.0 1 1 1 -999.0
21 28.0 529.000 680.000 400.0 5 5 5 3.3
22 20.0 27.660 115.000 148.0 5 5 5 11.0
23 3.9 0.120 1.000 16.0 3 1 2 -999.0
24 39.3 207.000 406.000 252.0 1 4 1 4.7
25 41.0 85.000 325.000 310.0 1 3 1 -999.0
26 16.2 36.330 119.500 63.0 1 1 1 10.4
27 9.0 0.101 4.000 28.0 5 1 3 7.4
28 7.6 1.040 5.500 68.0 5 3 4 2.1
29 46.0 521.000 655.000 336.0 5 5 5 -999.0
30 22.4 100.000 157.000 100.0 1 1 1 -999.0
31 16.3 35.000 56.000 33.0 3 5 4 7.7
32 2.6 0.005 0.140 21.5 5 2 4 17.9
33 24.0 0.010 0.250 50.0 1 1 1 6.1
34 100.0 62.000 1320.000 267.0 1 1 1 8.2
35 -999.0 0.122 3.000 30.0 2 1 1 8.4
36 -999.0 1.350 8.100 45.0 3 1 3 11.9
37 3.2 0.023 0.400 19.0 4 1 3 10.8
38 2.0 0.048 0.330 30.0 4 1 3 13.8
39 5.0 1.700 6.300 12.0 2 1 1 14.3
40 6.5 3.500 10.800 120.0 2 1 1 -999.0
41 23.6 250.000 490.000 440.0 5 5 5 15.2
42 12.0 0.480 15.500 140.0 2 2 2 10.0
43 20.2 10.000 115.000 170.0 4 4 4 11.9
44 13.0 1.620 11.400 17.0 2 1 2 6.5
45 27.0 192.000 180.000 115.0 4 4 4 7.5
46 18.0 2.500 12.100 31.0 5 5 5 -999.0
47 13.7 4.288 39.200 63.0 2 2 2 10.6
48 4.7 0.280 1.900 21.0 3 1 3 7.4
49 9.8 4.235 50.400 52.0 1 1 1 8.4
50 29.0 6.800 179.000 164.0 2 3 2 5.7
51 7.0 0.750 12.300 225.0 2 2 2 4.9
52 6.0 3.600 21.000 225.0 3 2 3 -999.0
53 17.0 14.830 98.200 150.0 5 5 5 3.2
54 20.0 55.500 175.000 151.0 5 5 5 -999.0
55 12.7 1.400 12.500 90.0 2 2 2 8.1
56 3.5 0.060 1.000 -999.0 3 1 1 11.0
57 4.5 0.900 2.600 38.0 2 1 2 -999.0
58 7.5 2.000 17.000 200.0 3 1 3 13.2
59 2.3 0.104 2.500 46.0 3 2 2 9.7
60 13.0 4.050 58.000 210.0 4 3 4 12.8
61 3.0 3.500 3.900 14.0 2 1 1 4.9
62 24.0 4.190 12.300 60.0 3 1 2 -999.0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) L wb wbr tg P
-59.38472 0.66987 -0.15338 0.11272 -0.67714 1.93900
S `D\r`
-2.07827 -0.04529
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-940.74 12.58 63.24 77.95 317.97
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -59.38472 76.26155 -0.779 0.43956
L 0.66987 0.20058 3.340 0.00153 **
wb -0.15338 0.09240 -1.660 0.10274
wbr 0.11272 0.09734 1.158 0.25193
tg -0.67714 0.26466 -2.559 0.01335 *
P 1.93900 34.87417 0.056 0.95587
S -2.07827 39.64169 -0.052 0.95838
`D\r` -0.04529 0.09202 -0.492 0.62461
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 271.8 on 54 degrees of freedom
Multiple R-squared: 0.2871, Adjusted R-squared: 0.1947
F-statistic: 3.106 on 7 and 54 DF, p-value: 0.00793
> 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
+ }
[,1] [,2] [,3]
[1,] 0.9885825 2.283503e-02 1.141752e-02
[2,] 0.9772781 4.544384e-02 2.272192e-02
[3,] 0.9903321 1.933570e-02 9.667851e-03
[4,] 0.9997814 4.372155e-04 2.186077e-04
[5,] 0.9994820 1.035986e-03 5.179930e-04
[6,] 0.9988317 2.336555e-03 1.168277e-03
[7,] 0.9974947 5.010615e-03 2.505308e-03
[8,] 0.9949739 1.005225e-02 5.026127e-03
[9,] 0.9915409 1.691821e-02 8.459107e-03
[10,] 0.9896695 2.066099e-02 1.033050e-02
[11,] 0.9864446 2.711088e-02 1.355544e-02
[12,] 0.9766600 4.668000e-02 2.334000e-02
[13,] 0.9634815 7.303691e-02 3.651846e-02
[14,] 0.9436012 1.127976e-01 5.639881e-02
[15,] 0.9150052 1.699896e-01 8.499481e-02
[16,] 0.8796686 2.406629e-01 1.203314e-01
[17,] 0.8356917 3.286166e-01 1.643083e-01
[18,] 0.7791208 4.417584e-01 2.208792e-01
[19,] 0.7327079 5.345842e-01 2.672921e-01
[20,] 0.6593896 6.812208e-01 3.406104e-01
[21,] 0.5834770 8.330461e-01 4.165230e-01
[22,] 0.5095824 9.808352e-01 4.904176e-01
[23,] 0.4320472 8.640944e-01 5.679528e-01
[24,] 0.3868188 7.736375e-01 6.131812e-01
[25,] 0.9821061 3.578786e-02 1.789393e-02
[26,] 1.0000000 1.324723e-28 6.623616e-29
[27,] 1.0000000 7.849893e-27 3.924947e-27
[28,] 1.0000000 4.546449e-25 2.273224e-25
[29,] 1.0000000 2.436430e-23 1.218215e-23
[30,] 1.0000000 1.173188e-21 5.865939e-22
[31,] 1.0000000 8.727041e-21 4.363520e-21
[32,] 1.0000000 4.682625e-19 2.341312e-19
[33,] 1.0000000 2.489244e-17 1.244622e-17
[34,] 1.0000000 8.734445e-16 4.367222e-16
[35,] 1.0000000 1.219173e-14 6.095865e-15
[36,] 1.0000000 4.997260e-13 2.498630e-13
[37,] 1.0000000 1.903925e-11 9.519625e-12
[38,] 1.0000000 1.004985e-09 5.024923e-10
[39,] 1.0000000 4.885821e-08 2.442910e-08
[40,] 0.9999992 1.623292e-06 8.116458e-07
[41,] 0.9999653 6.945045e-05 3.472523e-05
> postscript(file="/var/www/rcomp/tmp/13pt61292869839.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()
null device
1
> postscript(file="/var/www/rcomp/tmp/23pt61292869839.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()
null device
1
> postscript(file="/var/www/rcomp/tmp/3wgsr1292869839.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()
null device
1
> postscript(file="/var/www/rcomp/tmp/4wgsr1292869839.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()
null device
1
> postscript(file="/var/www/rcomp/tmp/5wgsr1292869839.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()
null device
1
> (myerror <- as.ts(mysum$resid))
Time Series:
Start = 1
End = 62
Frequency = 1
1 2 3 4 5
-147.30571345 86.91671712 -912.67742455 -259.39333369 314.34712725
6 7 8 9 10
154.40625354 86.31138852 317.96621337 94.07152003 148.52066164
11 12 13 14 15
139.67826344 257.32650138 57.82451592 -736.86485531 62.29009405
16 17 18 19 20
109.98502930 62.72442143 58.74124435 207.15385159 106.80121766
21 22 23 24 25
-203.53627196 66.39188596 20.56304831 -10.90057520 10.04799065
26 27 28 29 30
63.76283831 73.79098793 65.44344886 -222.18709032 -16.81774264
31 32 33 34 35
58.10683074 68.21071954 78.87294293 291.40083837 -940.74457255
36 37 38 39 40
-937.48359276 67.98240653 65.65084369 65.00068840 7.92354455
41 42 43 44 45
-54.15475272 59.74518805 75.16621732 74.99793454 -23.97975128
46 47 48 49 50
32.91351340 71.23607507 68.48376253 69.41329039 92.74961729
51 52 53 54 55
44.26124133 -0.02007406 68.83089258 10.86897717 64.91881118
56 57 58 59 60
178.27545564 17.72894266 52.53308979 59.56253551 60.68037471
61 62
60.77569866 34.71009731
> postscript(file="/var/www/rcomp/tmp/67q9c1292869839.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 62
Frequency = 1
lag(myerror, k = 1) myerror
0 -147.30571345 NA
1 86.91671712 -147.30571345
2 -912.67742455 86.91671712
3 -259.39333369 -912.67742455
4 314.34712725 -259.39333369
5 154.40625354 314.34712725
6 86.31138852 154.40625354
7 317.96621337 86.31138852
8 94.07152003 317.96621337
9 148.52066164 94.07152003
10 139.67826344 148.52066164
11 257.32650138 139.67826344
12 57.82451592 257.32650138
13 -736.86485531 57.82451592
14 62.29009405 -736.86485531
15 109.98502930 62.29009405
16 62.72442143 109.98502930
17 58.74124435 62.72442143
18 207.15385159 58.74124435
19 106.80121766 207.15385159
20 -203.53627196 106.80121766
21 66.39188596 -203.53627196
22 20.56304831 66.39188596
23 -10.90057520 20.56304831
24 10.04799065 -10.90057520
25 63.76283831 10.04799065
26 73.79098793 63.76283831
27 65.44344886 73.79098793
28 -222.18709032 65.44344886
29 -16.81774264 -222.18709032
30 58.10683074 -16.81774264
31 68.21071954 58.10683074
32 78.87294293 68.21071954
33 291.40083837 78.87294293
34 -940.74457255 291.40083837
35 -937.48359276 -940.74457255
36 67.98240653 -937.48359276
37 65.65084369 67.98240653
38 65.00068840 65.65084369
39 7.92354455 65.00068840
40 -54.15475272 7.92354455
41 59.74518805 -54.15475272
42 75.16621732 59.74518805
43 74.99793454 75.16621732
44 -23.97975128 74.99793454
45 32.91351340 -23.97975128
46 71.23607507 32.91351340
47 68.48376253 71.23607507
48 69.41329039 68.48376253
49 92.74961729 69.41329039
50 44.26124133 92.74961729
51 -0.02007406 44.26124133
52 68.83089258 -0.02007406
53 10.86897717 68.83089258
54 64.91881118 10.86897717
55 178.27545564 64.91881118
56 17.72894266 178.27545564
57 52.53308979 17.72894266
58 59.56253551 52.53308979
59 60.68037471 59.56253551
60 60.77569866 60.68037471
61 34.71009731 60.77569866
62 NA 34.71009731
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 86.91671712 -147.30571345
[2,] -912.67742455 86.91671712
[3,] -259.39333369 -912.67742455
[4,] 314.34712725 -259.39333369
[5,] 154.40625354 314.34712725
[6,] 86.31138852 154.40625354
[7,] 317.96621337 86.31138852
[8,] 94.07152003 317.96621337
[9,] 148.52066164 94.07152003
[10,] 139.67826344 148.52066164
[11,] 257.32650138 139.67826344
[12,] 57.82451592 257.32650138
[13,] -736.86485531 57.82451592
[14,] 62.29009405 -736.86485531
[15,] 109.98502930 62.29009405
[16,] 62.72442143 109.98502930
[17,] 58.74124435 62.72442143
[18,] 207.15385159 58.74124435
[19,] 106.80121766 207.15385159
[20,] -203.53627196 106.80121766
[21,] 66.39188596 -203.53627196
[22,] 20.56304831 66.39188596
[23,] -10.90057520 20.56304831
[24,] 10.04799065 -10.90057520
[25,] 63.76283831 10.04799065
[26,] 73.79098793 63.76283831
[27,] 65.44344886 73.79098793
[28,] -222.18709032 65.44344886
[29,] -16.81774264 -222.18709032
[30,] 58.10683074 -16.81774264
[31,] 68.21071954 58.10683074
[32,] 78.87294293 68.21071954
[33,] 291.40083837 78.87294293
[34,] -940.74457255 291.40083837
[35,] -937.48359276 -940.74457255
[36,] 67.98240653 -937.48359276
[37,] 65.65084369 67.98240653
[38,] 65.00068840 65.65084369
[39,] 7.92354455 65.00068840
[40,] -54.15475272 7.92354455
[41,] 59.74518805 -54.15475272
[42,] 75.16621732 59.74518805
[43,] 74.99793454 75.16621732
[44,] -23.97975128 74.99793454
[45,] 32.91351340 -23.97975128
[46,] 71.23607507 32.91351340
[47,] 68.48376253 71.23607507
[48,] 69.41329039 68.48376253
[49,] 92.74961729 69.41329039
[50,] 44.26124133 92.74961729
[51,] -0.02007406 44.26124133
[52,] 68.83089258 -0.02007406
[53,] 10.86897717 68.83089258
[54,] 64.91881118 10.86897717
[55,] 178.27545564 64.91881118
[56,] 17.72894266 178.27545564
[57,] 52.53308979 17.72894266
[58,] 59.56253551 52.53308979
[59,] 60.68037471 59.56253551
[60,] 60.77569866 60.68037471
[61,] 34.71009731 60.77569866
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 86.91671712 -147.30571345
2 -912.67742455 86.91671712
3 -259.39333369 -912.67742455
4 314.34712725 -259.39333369
5 154.40625354 314.34712725
6 86.31138852 154.40625354
7 317.96621337 86.31138852
8 94.07152003 317.96621337
9 148.52066164 94.07152003
10 139.67826344 148.52066164
11 257.32650138 139.67826344
12 57.82451592 257.32650138
13 -736.86485531 57.82451592
14 62.29009405 -736.86485531
15 109.98502930 62.29009405
16 62.72442143 109.98502930
17 58.74124435 62.72442143
18 207.15385159 58.74124435
19 106.80121766 207.15385159
20 -203.53627196 106.80121766
21 66.39188596 -203.53627196
22 20.56304831 66.39188596
23 -10.90057520 20.56304831
24 10.04799065 -10.90057520
25 63.76283831 10.04799065
26 73.79098793 63.76283831
27 65.44344886 73.79098793
28 -222.18709032 65.44344886
29 -16.81774264 -222.18709032
30 58.10683074 -16.81774264
31 68.21071954 58.10683074
32 78.87294293 68.21071954
33 291.40083837 78.87294293
34 -940.74457255 291.40083837
35 -937.48359276 -940.74457255
36 67.98240653 -937.48359276
37 65.65084369 67.98240653
38 65.00068840 65.65084369
39 7.92354455 65.00068840
40 -54.15475272 7.92354455
41 59.74518805 -54.15475272
42 75.16621732 59.74518805
43 74.99793454 75.16621732
44 -23.97975128 74.99793454
45 32.91351340 -23.97975128
46 71.23607507 32.91351340
47 68.48376253 71.23607507
48 69.41329039 68.48376253
49 92.74961729 69.41329039
50 44.26124133 92.74961729
51 -0.02007406 44.26124133
52 68.83089258 -0.02007406
53 10.86897717 68.83089258
54 64.91881118 10.86897717
55 178.27545564 64.91881118
56 17.72894266 178.27545564
57 52.53308979 17.72894266
58 59.56253551 52.53308979
59 60.68037471 59.56253551
60 60.77569866 60.68037471
61 34.71009731 60.77569866
> 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()
null device
1
> postscript(file="/var/www/rcomp/tmp/77q9c1292869839.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()
null device
1
> postscript(file="/var/www/rcomp/tmp/8hz9x1292869839.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()
null device
1
> postscript(file="/var/www/rcomp/tmp/9hz9x1292869839.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()
null device
1
> if (n > n25) {
+ postscript(file="/var/www/rcomp/tmp/10s88z1292869839.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()
+ }
null device
1
>
> #Note: the /var/www/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/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/www/rcomp/tmp/11vqon1292869839.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
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/www/rcomp/tmp/12oi5q1292869839.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/www/rcomp/tmp/13d1321292869839.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
Forecast', 1, TRUE)
> a<-table.element(a, 'Residuals
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/www/rcomp/tmp/14ns251292869839.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/www/rcomp/tmp/15rtib1292869839.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/www/rcomp/tmp/16nkyj1292869839.tab")
+ }
>
> try(system("convert tmp/13pt61292869839.ps tmp/13pt61292869839.png",intern=TRUE))
character(0)
> try(system("convert tmp/23pt61292869839.ps tmp/23pt61292869839.png",intern=TRUE))
character(0)
> try(system("convert tmp/3wgsr1292869839.ps tmp/3wgsr1292869839.png",intern=TRUE))
character(0)
> try(system("convert tmp/4wgsr1292869839.ps tmp/4wgsr1292869839.png",intern=TRUE))
character(0)
> try(system("convert tmp/5wgsr1292869839.ps tmp/5wgsr1292869839.png",intern=TRUE))
character(0)
> try(system("convert tmp/67q9c1292869839.ps tmp/67q9c1292869839.png",intern=TRUE))
character(0)
> try(system("convert tmp/77q9c1292869839.ps tmp/77q9c1292869839.png",intern=TRUE))
character(0)
> try(system("convert tmp/8hz9x1292869839.ps tmp/8hz9x1292869839.png",intern=TRUE))
character(0)
> try(system("convert tmp/9hz9x1292869839.ps tmp/9hz9x1292869839.png",intern=TRUE))
character(0)
> try(system("convert tmp/10s88z1292869839.ps tmp/10s88z1292869839.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.170 0.860 4.018