R version 2.6.2 (2008-02-08)
Copyright (C) 2008 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
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.
Natural language support but running in an English locale
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(56421,53536,52408,41454,38271,35306,26414,31917,38030,27534,18387,50556,43901,43899,37532,40357,35489,29027,34485,42598,30306,26451,47460,50104,61465,53726,39477,43895,31481,29896,33842,39120,33702,25094,51442,45594,52518,48564,41745,49585,32747,33379,35645,37034,35681,20972,58552,54955,51570,51145,46641,35704,33253,35193,41668,34865,21210,56126,49231,59723,48103,47472,50497,40059,34149,36860,46356,36577),dim=c(1,68),dimnames=list(c(''),1:68))
> y <- array(NA,dim=c(1,68),dimnames=list(c(''),1:68))
> 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'
> #'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
M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 56421 1 0 0 0 0 0 0 0 0 0 0 1
2 53536 0 1 0 0 0 0 0 0 0 0 0 2
3 52408 0 0 1 0 0 0 0 0 0 0 0 3
4 41454 0 0 0 1 0 0 0 0 0 0 0 4
5 38271 0 0 0 0 1 0 0 0 0 0 0 5
6 35306 0 0 0 0 0 1 0 0 0 0 0 6
7 26414 0 0 0 0 0 0 1 0 0 0 0 7
8 31917 0 0 0 0 0 0 0 1 0 0 0 8
9 38030 0 0 0 0 0 0 0 0 1 0 0 9
10 27534 0 0 0 0 0 0 0 0 0 1 0 10
11 18387 0 0 0 0 0 0 0 0 0 0 1 11
12 50556 0 0 0 0 0 0 0 0 0 0 0 12
13 43901 1 0 0 0 0 0 0 0 0 0 0 13
14 43899 0 1 0 0 0 0 0 0 0 0 0 14
15 37532 0 0 1 0 0 0 0 0 0 0 0 15
16 40357 0 0 0 1 0 0 0 0 0 0 0 16
17 35489 0 0 0 0 1 0 0 0 0 0 0 17
18 29027 0 0 0 0 0 1 0 0 0 0 0 18
19 34485 0 0 0 0 0 0 1 0 0 0 0 19
20 42598 0 0 0 0 0 0 0 1 0 0 0 20
21 30306 0 0 0 0 0 0 0 0 1 0 0 21
22 26451 0 0 0 0 0 0 0 0 0 1 0 22
23 47460 0 0 0 0 0 0 0 0 0 0 1 23
24 50104 0 0 0 0 0 0 0 0 0 0 0 24
25 61465 1 0 0 0 0 0 0 0 0 0 0 25
26 53726 0 1 0 0 0 0 0 0 0 0 0 26
27 39477 0 0 1 0 0 0 0 0 0 0 0 27
28 43895 0 0 0 1 0 0 0 0 0 0 0 28
29 31481 0 0 0 0 1 0 0 0 0 0 0 29
30 29896 0 0 0 0 0 1 0 0 0 0 0 30
31 33842 0 0 0 0 0 0 1 0 0 0 0 31
32 39120 0 0 0 0 0 0 0 1 0 0 0 32
33 33702 0 0 0 0 0 0 0 0 1 0 0 33
34 25094 0 0 0 0 0 0 0 0 0 1 0 34
35 51442 0 0 0 0 0 0 0 0 0 0 1 35
36 45594 0 0 0 0 0 0 0 0 0 0 0 36
37 52518 1 0 0 0 0 0 0 0 0 0 0 37
38 48564 0 1 0 0 0 0 0 0 0 0 0 38
39 41745 0 0 1 0 0 0 0 0 0 0 0 39
40 49585 0 0 0 1 0 0 0 0 0 0 0 40
41 32747 0 0 0 0 1 0 0 0 0 0 0 41
42 33379 0 0 0 0 0 1 0 0 0 0 0 42
43 35645 0 0 0 0 0 0 1 0 0 0 0 43
44 37034 0 0 0 0 0 0 0 1 0 0 0 44
45 35681 0 0 0 0 0 0 0 0 1 0 0 45
46 20972 0 0 0 0 0 0 0 0 0 1 0 46
47 58552 0 0 0 0 0 0 0 0 0 0 1 47
48 54955 0 0 0 0 0 0 0 0 0 0 0 48
49 51570 1 0 0 0 0 0 0 0 0 0 0 49
50 51145 0 1 0 0 0 0 0 0 0 0 0 50
51 46641 0 0 1 0 0 0 0 0 0 0 0 51
52 35704 0 0 0 1 0 0 0 0 0 0 0 52
53 33253 0 0 0 0 1 0 0 0 0 0 0 53
54 35193 0 0 0 0 0 1 0 0 0 0 0 54
55 41668 0 0 0 0 0 0 1 0 0 0 0 55
56 34865 0 0 0 0 0 0 0 1 0 0 0 56
57 21210 0 0 0 0 0 0 0 0 1 0 0 57
58 56126 0 0 0 0 0 0 0 0 0 1 0 58
59 49231 0 0 0 0 0 0 0 0 0 0 1 59
60 59723 0 0 0 0 0 0 0 0 0 0 0 60
61 48103 1 0 0 0 0 0 0 0 0 0 0 61
62 47472 0 1 0 0 0 0 0 0 0 0 0 62
63 50497 0 0 1 0 0 0 0 0 0 0 0 63
64 40059 0 0 0 1 0 0 0 0 0 0 0 64
65 34149 0 0 0 0 1 0 0 0 0 0 0 65
66 36860 0 0 0 0 0 1 0 0 0 0 0 66
67 46356 0 0 0 0 0 0 1 0 0 0 0 67
68 36577 0 0 0 0 0 0 0 1 0 0 0 68
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) M1 M2 M3 M4 M5
49460.50 521.86 -2159.86 -7242.57 -10192.63 -17879.01
M6 M7 M8 M9 M10 M11
-18909.57 -15860.45 -15319.34 -20173.44 -20799.56 -7096.28
t
75.72
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-24810.1 -3051.4 -395.5 3392.9 23073.3
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 49460.50 3587.51 13.787 < 2e-16 ***
M1 521.86 4350.40 0.120 0.904954
M2 -2159.86 4348.35 -0.497 0.621375
M3 -7242.57 4346.75 -1.666 0.101359
M4 -10192.63 4345.61 -2.346 0.022636 *
M5 -17879.01 4344.92 -4.115 0.000131 ***
M6 -18909.57 4344.69 -4.352 5.90e-05 ***
M7 -15860.45 4344.92 -3.650 0.000584 ***
M8 -15319.34 4345.61 -3.525 0.000861 ***
M9 -20173.44 4539.85 -4.444 4.33e-05 ***
M10 -20799.56 4538.76 -4.583 2.68e-05 ***
M11 -7096.28 4538.10 -1.564 0.123621
t 75.72 44.57 1.699 0.094963 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 7175 on 55 degrees of freedom
Multiple R-squared: 0.5808, Adjusted R-squared: 0.4894
F-statistic: 6.351 on 12 and 55 DF, p-value: 6.806e-07
> 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.2882178 0.576435648 0.7117821760
[2,] 0.2315709 0.463141706 0.7684291468
[3,] 0.1292436 0.258487111 0.8707564447
[4,] 0.3706799 0.741359742 0.6293201291
[5,] 0.5704681 0.859063852 0.4295319262
[6,] 0.4715272 0.943054451 0.5284727743
[7,] 0.3789959 0.757991893 0.6210040533
[8,] 0.9001428 0.199714491 0.0998572453
[9,] 0.8499404 0.300119234 0.1500596172
[10,] 0.8846452 0.230709641 0.1153548206
[11,] 0.8531813 0.293637450 0.1468187248
[12,] 0.8198481 0.360303814 0.1801519069
[13,] 0.7655479 0.468904118 0.2344520591
[14,] 0.7119441 0.576111730 0.2880558649
[15,] 0.6363539 0.727292296 0.3636461478
[16,] 0.5662890 0.867421930 0.4337109651
[17,] 0.4979024 0.995804792 0.5020976038
[18,] 0.4427217 0.885443423 0.5572782883
[19,] 0.4224337 0.844867482 0.5775662588
[20,] 0.5076432 0.984713548 0.4923567738
[21,] 0.4955085 0.991017084 0.5044914580
[22,] 0.4221060 0.844211933 0.5778940335
[23,] 0.3415195 0.683039027 0.6584804865
[24,] 0.2788128 0.557625520 0.7211872400
[25,] 0.3311019 0.662203783 0.6688981086
[26,] 0.2594657 0.518931354 0.7405343230
[27,] 0.1899315 0.379862960 0.8100685199
[28,] 0.1463355 0.292671072 0.8536644639
[29,] 0.1039509 0.207901716 0.8960491422
[30,] 0.1441642 0.288328321 0.8558358396
[31,] 0.9769927 0.046014560 0.0230072798
[32,] 0.9981270 0.003746074 0.0018730369
[33,] 0.9963465 0.007306938 0.0036534692
[34,] 0.9964499 0.007100143 0.0035500717
[35,] 0.9994782 0.001043653 0.0005218265
[36,] 0.9973866 0.005226894 0.0026134472
[37,] 0.9918030 0.016393908 0.0081969542
> postscript(file="/var/www/html/rcomp/tmp/1ggbt1210263175.ps",horizontal=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/html/rcomp/tmp/2gegb1210263175.ps",horizontal=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/html/rcomp/tmp/3onsg1210263175.ps",horizontal=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/html/rcomp/tmp/4kw8o1210263175.ps",horizontal=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/html/rcomp/tmp/5ry781210263175.ps",horizontal=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 = 68
Frequency = 1
1 2 3 4 5 6
6362.91667 6083.91667 9962.91667 1883.25000 6310.91667 4300.75000
7 8 9 10 11 12
-7716.08333 -2829.91667 8061.46667 -1884.13333 -24810.13333 186.86667
13 14 15 16 17 18
-7065.71667 -4461.71667 -5821.71667 -122.38333 2620.28333 -2886.88333
19 20 21 22 23 24
-553.71667 6942.45000 -571.16667 -3875.76667 3354.23333 -1173.76667
25 26 27 28 29 30
9589.65000 4456.65000 -4785.35000 2506.98333 -2296.35000 -2926.51667
31 32 33 34 35 36
-2105.35000 2555.81667 1916.20000 -6141.40000 6427.60000 -6592.40000
37 38 39 40 41 42
-265.98333 -1613.98333 -3425.98333 7288.35000 -1938.98333 -352.15000
43 44 45 46 47 48
-1210.98333 -438.81667 2986.56667 -11172.03333 12628.96667 1859.96667
49 50 51 52 53 54
-2122.61667 58.38333 561.38333 -7501.28333 -2341.61667 553.21667
55 56 57 58 59 60
3903.38333 -3516.45000 -12393.06667 23073.33333 2399.33333 5719.33333
61 62 63 64 65 66
-6498.25000 -4523.25000 3508.75000 -4054.91667 -2354.25000 1311.58333
67 68
7682.75000 -2713.08333
> postscript(file="/var/www/html/rcomp/tmp/6ffyp1210263175.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 68
Frequency = 1
lag(myerror, k = 1) myerror
0 6362.91667 NA
1 6083.91667 6362.91667
2 9962.91667 6083.91667
3 1883.25000 9962.91667
4 6310.91667 1883.25000
5 4300.75000 6310.91667
6 -7716.08333 4300.75000
7 -2829.91667 -7716.08333
8 8061.46667 -2829.91667
9 -1884.13333 8061.46667
10 -24810.13333 -1884.13333
11 186.86667 -24810.13333
12 -7065.71667 186.86667
13 -4461.71667 -7065.71667
14 -5821.71667 -4461.71667
15 -122.38333 -5821.71667
16 2620.28333 -122.38333
17 -2886.88333 2620.28333
18 -553.71667 -2886.88333
19 6942.45000 -553.71667
20 -571.16667 6942.45000
21 -3875.76667 -571.16667
22 3354.23333 -3875.76667
23 -1173.76667 3354.23333
24 9589.65000 -1173.76667
25 4456.65000 9589.65000
26 -4785.35000 4456.65000
27 2506.98333 -4785.35000
28 -2296.35000 2506.98333
29 -2926.51667 -2296.35000
30 -2105.35000 -2926.51667
31 2555.81667 -2105.35000
32 1916.20000 2555.81667
33 -6141.40000 1916.20000
34 6427.60000 -6141.40000
35 -6592.40000 6427.60000
36 -265.98333 -6592.40000
37 -1613.98333 -265.98333
38 -3425.98333 -1613.98333
39 7288.35000 -3425.98333
40 -1938.98333 7288.35000
41 -352.15000 -1938.98333
42 -1210.98333 -352.15000
43 -438.81667 -1210.98333
44 2986.56667 -438.81667
45 -11172.03333 2986.56667
46 12628.96667 -11172.03333
47 1859.96667 12628.96667
48 -2122.61667 1859.96667
49 58.38333 -2122.61667
50 561.38333 58.38333
51 -7501.28333 561.38333
52 -2341.61667 -7501.28333
53 553.21667 -2341.61667
54 3903.38333 553.21667
55 -3516.45000 3903.38333
56 -12393.06667 -3516.45000
57 23073.33333 -12393.06667
58 2399.33333 23073.33333
59 5719.33333 2399.33333
60 -6498.25000 5719.33333
61 -4523.25000 -6498.25000
62 3508.75000 -4523.25000
63 -4054.91667 3508.75000
64 -2354.25000 -4054.91667
65 1311.58333 -2354.25000
66 7682.75000 1311.58333
67 -2713.08333 7682.75000
68 NA -2713.08333
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 6083.91667 6362.91667
[2,] 9962.91667 6083.91667
[3,] 1883.25000 9962.91667
[4,] 6310.91667 1883.25000
[5,] 4300.75000 6310.91667
[6,] -7716.08333 4300.75000
[7,] -2829.91667 -7716.08333
[8,] 8061.46667 -2829.91667
[9,] -1884.13333 8061.46667
[10,] -24810.13333 -1884.13333
[11,] 186.86667 -24810.13333
[12,] -7065.71667 186.86667
[13,] -4461.71667 -7065.71667
[14,] -5821.71667 -4461.71667
[15,] -122.38333 -5821.71667
[16,] 2620.28333 -122.38333
[17,] -2886.88333 2620.28333
[18,] -553.71667 -2886.88333
[19,] 6942.45000 -553.71667
[20,] -571.16667 6942.45000
[21,] -3875.76667 -571.16667
[22,] 3354.23333 -3875.76667
[23,] -1173.76667 3354.23333
[24,] 9589.65000 -1173.76667
[25,] 4456.65000 9589.65000
[26,] -4785.35000 4456.65000
[27,] 2506.98333 -4785.35000
[28,] -2296.35000 2506.98333
[29,] -2926.51667 -2296.35000
[30,] -2105.35000 -2926.51667
[31,] 2555.81667 -2105.35000
[32,] 1916.20000 2555.81667
[33,] -6141.40000 1916.20000
[34,] 6427.60000 -6141.40000
[35,] -6592.40000 6427.60000
[36,] -265.98333 -6592.40000
[37,] -1613.98333 -265.98333
[38,] -3425.98333 -1613.98333
[39,] 7288.35000 -3425.98333
[40,] -1938.98333 7288.35000
[41,] -352.15000 -1938.98333
[42,] -1210.98333 -352.15000
[43,] -438.81667 -1210.98333
[44,] 2986.56667 -438.81667
[45,] -11172.03333 2986.56667
[46,] 12628.96667 -11172.03333
[47,] 1859.96667 12628.96667
[48,] -2122.61667 1859.96667
[49,] 58.38333 -2122.61667
[50,] 561.38333 58.38333
[51,] -7501.28333 561.38333
[52,] -2341.61667 -7501.28333
[53,] 553.21667 -2341.61667
[54,] 3903.38333 553.21667
[55,] -3516.45000 3903.38333
[56,] -12393.06667 -3516.45000
[57,] 23073.33333 -12393.06667
[58,] 2399.33333 23073.33333
[59,] 5719.33333 2399.33333
[60,] -6498.25000 5719.33333
[61,] -4523.25000 -6498.25000
[62,] 3508.75000 -4523.25000
[63,] -4054.91667 3508.75000
[64,] -2354.25000 -4054.91667
[65,] 1311.58333 -2354.25000
[66,] 7682.75000 1311.58333
[67,] -2713.08333 7682.75000
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 6083.91667 6362.91667
2 9962.91667 6083.91667
3 1883.25000 9962.91667
4 6310.91667 1883.25000
5 4300.75000 6310.91667
6 -7716.08333 4300.75000
7 -2829.91667 -7716.08333
8 8061.46667 -2829.91667
9 -1884.13333 8061.46667
10 -24810.13333 -1884.13333
11 186.86667 -24810.13333
12 -7065.71667 186.86667
13 -4461.71667 -7065.71667
14 -5821.71667 -4461.71667
15 -122.38333 -5821.71667
16 2620.28333 -122.38333
17 -2886.88333 2620.28333
18 -553.71667 -2886.88333
19 6942.45000 -553.71667
20 -571.16667 6942.45000
21 -3875.76667 -571.16667
22 3354.23333 -3875.76667
23 -1173.76667 3354.23333
24 9589.65000 -1173.76667
25 4456.65000 9589.65000
26 -4785.35000 4456.65000
27 2506.98333 -4785.35000
28 -2296.35000 2506.98333
29 -2926.51667 -2296.35000
30 -2105.35000 -2926.51667
31 2555.81667 -2105.35000
32 1916.20000 2555.81667
33 -6141.40000 1916.20000
34 6427.60000 -6141.40000
35 -6592.40000 6427.60000
36 -265.98333 -6592.40000
37 -1613.98333 -265.98333
38 -3425.98333 -1613.98333
39 7288.35000 -3425.98333
40 -1938.98333 7288.35000
41 -352.15000 -1938.98333
42 -1210.98333 -352.15000
43 -438.81667 -1210.98333
44 2986.56667 -438.81667
45 -11172.03333 2986.56667
46 12628.96667 -11172.03333
47 1859.96667 12628.96667
48 -2122.61667 1859.96667
49 58.38333 -2122.61667
50 561.38333 58.38333
51 -7501.28333 561.38333
52 -2341.61667 -7501.28333
53 553.21667 -2341.61667
54 3903.38333 553.21667
55 -3516.45000 3903.38333
56 -12393.06667 -3516.45000
57 23073.33333 -12393.06667
58 2399.33333 23073.33333
59 5719.33333 2399.33333
60 -6498.25000 5719.33333
61 -4523.25000 -6498.25000
62 3508.75000 -4523.25000
63 -4054.91667 3508.75000
64 -2354.25000 -4054.91667
65 1311.58333 -2354.25000
66 7682.75000 1311.58333
67 -2713.08333 7682.75000
> 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/html/rcomp/tmp/7yrx91210263175.ps",horizontal=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/html/rcomp/tmp/8biqf1210263175.ps",horizontal=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/html/rcomp/tmp/9acnj1210263175.ps",horizontal=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/html/rcomp/tmp/107bxm1210263175.ps",horizontal=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/html/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/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/html/rcomp/tmp/11si5h1210263176.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/html/rcomp/tmp/12nfgt1210263176.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/html/rcomp/tmp/13nz8e1210263176.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/html/rcomp/tmp/14h2uw1210263176.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/html/rcomp/tmp/15m28h1210263176.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/html/rcomp/tmp/16jqz11210263176.tab")
+ }
>
> system("convert tmp/1ggbt1210263175.ps tmp/1ggbt1210263175.png")
> system("convert tmp/2gegb1210263175.ps tmp/2gegb1210263175.png")
> system("convert tmp/3onsg1210263175.ps tmp/3onsg1210263175.png")
> system("convert tmp/4kw8o1210263175.ps tmp/4kw8o1210263175.png")
> system("convert tmp/5ry781210263175.ps tmp/5ry781210263175.png")
> system("convert tmp/6ffyp1210263175.ps tmp/6ffyp1210263175.png")
> system("convert tmp/7yrx91210263175.ps tmp/7yrx91210263175.png")
> system("convert tmp/8biqf1210263175.ps tmp/8biqf1210263175.png")
> system("convert tmp/9acnj1210263175.ps tmp/9acnj1210263175.png")
> system("convert tmp/107bxm1210263175.ps tmp/107bxm1210263175.png")
>
>
> proc.time()
user system elapsed
5.351 2.780 5.713