R version 2.9.0 (2009-04-17)
Copyright (C) 2009 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.
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(10519.20
+ ,1154.80
+ ,10414.90
+ ,1206.70
+ ,12476.80
+ ,1199.00
+ ,12384.60
+ ,1265.00
+ ,12266.70
+ ,1247.10
+ ,12919.90
+ ,1116.50
+ ,11497.30
+ ,1153.90
+ ,12142.00
+ ,1077.40
+ ,13919.40
+ ,1132.50
+ ,12656.80
+ ,1058.80
+ ,12034.10
+ ,1195.10
+ ,13199.70
+ ,1263.40
+ ,10881.30
+ ,1023.10
+ ,11301.20
+ ,1141.00
+ ,13643.90
+ ,1116.30
+ ,12517.00
+ ,1135.60
+ ,13981.10
+ ,1210.50
+ ,14275.70
+ ,1230.00
+ ,13425.00
+ ,1136.50
+ ,13565.70
+ ,1068.70
+ ,16216.30
+ ,1372.50
+ ,12970.00
+ ,1049.90
+ ,14079.90
+ ,1302.20
+ ,14235.00
+ ,1305.90
+ ,12213.40
+ ,1173.50
+ ,12581.00
+ ,1277.40
+ ,14130.40
+ ,1238.60
+ ,14210.80
+ ,1508.60
+ ,14378.50
+ ,1423.40
+ ,13142.80
+ ,1375.10
+ ,13714.70
+ ,1344.10
+ ,13621.90
+ ,1287.50
+ ,15379.80
+ ,1446.90
+ ,13306.30
+ ,1451.00
+ ,14391.20
+ ,1604.40
+ ,14909.90
+ ,1501.50
+ ,14025.40
+ ,1522.80
+ ,12951.20
+ ,1328.00
+ ,14344.30
+ ,1420.50
+ ,16093.40
+ ,1648.00
+ ,15413.60
+ ,1631.10
+ ,14705.70
+ ,1396.60
+ ,15972.80
+ ,1663.40
+ ,16241.40
+ ,1283.00
+ ,16626.40
+ ,1582.40
+ ,17136.20
+ ,1785.20
+ ,15622.90
+ ,1853.60
+ ,18003.90
+ ,1994.10
+ ,16136.10
+ ,2042.80
+ ,14423.70
+ ,1586.10
+ ,16789.40
+ ,1942.40
+ ,16782.20
+ ,1763.60
+ ,14133.80
+ ,1819.90
+ ,12607.00
+ ,1836.00
+ ,12004.50
+ ,1447.50
+ ,12175.40
+ ,1509.50
+ ,13268.00
+ ,1661.20
+ ,12299.30
+ ,1456.20
+ ,11800.60
+ ,1310.90
+ ,13873.30
+ ,1542.10
+ ,12315.00
+ ,1537.70)
+ ,dim=c(2
+ ,61)
+ ,dimnames=list(c('InvoerEU'
+ ,'InvoerAM')
+ ,1:61))
> y <- array(NA,dim=c(2,61),dimnames=list(c('InvoerEU','InvoerAM'),1:61))
> 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
Attaching package: 'zoo'
The following object(s) are masked from package:base :
as.Date.numeric
> 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
InvoerEU InvoerAM M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 10519.2 1154.8 1 0 0 0 0 0 0 0 0 0 0 1
2 10414.9 1206.7 0 1 0 0 0 0 0 0 0 0 0 2
3 12476.8 1199.0 0 0 1 0 0 0 0 0 0 0 0 3
4 12384.6 1265.0 0 0 0 1 0 0 0 0 0 0 0 4
5 12266.7 1247.1 0 0 0 0 1 0 0 0 0 0 0 5
6 12919.9 1116.5 0 0 0 0 0 1 0 0 0 0 0 6
7 11497.3 1153.9 0 0 0 0 0 0 1 0 0 0 0 7
8 12142.0 1077.4 0 0 0 0 0 0 0 1 0 0 0 8
9 13919.4 1132.5 0 0 0 0 0 0 0 0 1 0 0 9
10 12656.8 1058.8 0 0 0 0 0 0 0 0 0 1 0 10
11 12034.1 1195.1 0 0 0 0 0 0 0 0 0 0 1 11
12 13199.7 1263.4 0 0 0 0 0 0 0 0 0 0 0 12
13 10881.3 1023.1 1 0 0 0 0 0 0 0 0 0 0 13
14 11301.2 1141.0 0 1 0 0 0 0 0 0 0 0 0 14
15 13643.9 1116.3 0 0 1 0 0 0 0 0 0 0 0 15
16 12517.0 1135.6 0 0 0 1 0 0 0 0 0 0 0 16
17 13981.1 1210.5 0 0 0 0 1 0 0 0 0 0 0 17
18 14275.7 1230.0 0 0 0 0 0 1 0 0 0 0 0 18
19 13425.0 1136.5 0 0 0 0 0 0 1 0 0 0 0 19
20 13565.7 1068.7 0 0 0 0 0 0 0 1 0 0 0 20
21 16216.3 1372.5 0 0 0 0 0 0 0 0 1 0 0 21
22 12970.0 1049.9 0 0 0 0 0 0 0 0 0 1 0 22
23 14079.9 1302.2 0 0 0 0 0 0 0 0 0 0 1 23
24 14235.0 1305.9 0 0 0 0 0 0 0 0 0 0 0 24
25 12213.4 1173.5 1 0 0 0 0 0 0 0 0 0 0 25
26 12581.0 1277.4 0 1 0 0 0 0 0 0 0 0 0 26
27 14130.4 1238.6 0 0 1 0 0 0 0 0 0 0 0 27
28 14210.8 1508.6 0 0 0 1 0 0 0 0 0 0 0 28
29 14378.5 1423.4 0 0 0 0 1 0 0 0 0 0 0 29
30 13142.8 1375.1 0 0 0 0 0 1 0 0 0 0 0 30
31 13714.7 1344.1 0 0 0 0 0 0 1 0 0 0 0 31
32 13621.9 1287.5 0 0 0 0 0 0 0 1 0 0 0 32
33 15379.8 1446.9 0 0 0 0 0 0 0 0 1 0 0 33
34 13306.3 1451.0 0 0 0 0 0 0 0 0 0 1 0 34
35 14391.2 1604.4 0 0 0 0 0 0 0 0 0 0 1 35
36 14909.9 1501.5 0 0 0 0 0 0 0 0 0 0 0 36
37 14025.4 1522.8 1 0 0 0 0 0 0 0 0 0 0 37
38 12951.2 1328.0 0 1 0 0 0 0 0 0 0 0 0 38
39 14344.3 1420.5 0 0 1 0 0 0 0 0 0 0 0 39
40 16093.4 1648.0 0 0 0 1 0 0 0 0 0 0 0 40
41 15413.6 1631.1 0 0 0 0 1 0 0 0 0 0 0 41
42 14705.7 1396.6 0 0 0 0 0 1 0 0 0 0 0 42
43 15972.8 1663.4 0 0 0 0 0 0 1 0 0 0 0 43
44 16241.4 1283.0 0 0 0 0 0 0 0 1 0 0 0 44
45 16626.4 1582.4 0 0 0 0 0 0 0 0 1 0 0 45
46 17136.2 1785.2 0 0 0 0 0 0 0 0 0 1 0 46
47 15622.9 1853.6 0 0 0 0 0 0 0 0 0 0 1 47
48 18003.9 1994.1 0 0 0 0 0 0 0 0 0 0 0 48
49 16136.1 2042.8 1 0 0 0 0 0 0 0 0 0 0 49
50 14423.7 1586.1 0 1 0 0 0 0 0 0 0 0 0 50
51 16789.4 1942.4 0 0 1 0 0 0 0 0 0 0 0 51
52 16782.2 1763.6 0 0 0 1 0 0 0 0 0 0 0 52
53 14133.8 1819.9 0 0 0 0 1 0 0 0 0 0 0 53
54 12607.0 1836.0 0 0 0 0 0 1 0 0 0 0 0 54
55 12004.5 1447.5 0 0 0 0 0 0 1 0 0 0 0 55
56 12175.4 1509.5 0 0 0 0 0 0 0 1 0 0 0 56
57 13268.0 1661.2 0 0 0 0 0 0 0 0 1 0 0 57
58 12299.3 1456.2 0 0 0 0 0 0 0 0 0 1 0 58
59 11800.6 1310.9 0 0 0 0 0 0 0 0 0 0 1 59
60 13873.3 1542.1 0 0 0 0 0 0 0 0 0 0 0 60
61 12315.0 1537.7 1 0 0 0 0 0 0 0 0 0 0 61
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) InvoerAM M1 M2 M3 M4
7623.860 5.092 -1664.044 -1568.845 3.777 -272.409
M5 M6 M7 M8 M9 M10
-632.034 -737.127 -717.184 52.779 612.800 -379.126
M11 t
-926.162 -14.642
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-2838.7 -708.9 166.9 739.6 2675.4
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7623.860 1363.541 5.591 1.11e-06 ***
InvoerAM 5.092 1.075 4.735 2.05e-05 ***
M1 -1664.044 749.861 -2.219 0.0313 *
M2 -1568.845 792.083 -1.981 0.0535 .
M3 3.777 784.512 0.005 0.9962
M4 -272.409 783.344 -0.348 0.7296
M5 -632.034 782.186 -0.808 0.4231
M6 -737.127 784.080 -0.940 0.3520
M7 -717.184 790.500 -0.907 0.3689
M8 52.779 818.803 0.064 0.9489
M9 612.800 781.432 0.784 0.4369
M10 -379.126 793.576 -0.478 0.6350
M11 -926.162 781.651 -1.185 0.2420
t -14.642 15.071 -0.972 0.3362
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1232 on 47 degrees of freedom
Multiple R-squared: 0.5931, Adjusted R-squared: 0.4805
F-statistic: 5.269 on 13 and 47 DF, p-value: 1.079e-05
> 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,] 5.271206e-03 1.054241e-02 0.9947288
[2,] 1.707184e-02 3.414369e-02 0.9829282
[3,] 1.217708e-02 2.435417e-02 0.9878229
[4,] 3.606263e-03 7.212526e-03 0.9963937
[5,] 1.043814e-03 2.087629e-03 0.9989562
[6,] 1.183358e-03 2.366716e-03 0.9988166
[7,] 3.948618e-04 7.897236e-04 0.9996051
[8,] 1.449307e-04 2.898614e-04 0.9998551
[9,] 8.887492e-05 1.777498e-04 0.9999111
[10,] 3.702254e-05 7.404508e-05 0.9999630
[11,] 3.651125e-05 7.302250e-05 0.9999635
[12,] 4.600772e-05 9.201544e-05 0.9999540
[13,] 2.752032e-05 5.504063e-05 0.9999725
[14,] 1.031974e-03 2.063949e-03 0.9989680
[15,] 4.506531e-04 9.013062e-04 0.9995493
[16,] 3.558820e-04 7.117640e-04 0.9996441
[17,] 2.849317e-04 5.698634e-04 0.9997151
[18,] 4.105877e-04 8.211753e-04 0.9995894
[19,] 2.915112e-04 5.830223e-04 0.9997085
[20,] 3.333740e-04 6.667480e-04 0.9996666
[21,] 7.939142e-04 1.587828e-03 0.9992061
[22,] 1.411557e-03 2.823113e-03 0.9985884
[23,] 1.231261e-02 2.462522e-02 0.9876874
[24,] 4.034108e-01 8.068216e-01 0.5965892
[25,] 7.423139e-01 5.153722e-01 0.2576861
[26,] 8.957899e-01 2.084201e-01 0.1042101
[27,] 8.195317e-01 3.609365e-01 0.1804683
[28,] 7.414977e-01 5.170045e-01 0.2585023
> postscript(file="/var/www/html/rcomp/tmp/19cz91262182684.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/26d5x1262182684.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/30i061262182684.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/4qb2a1262182684.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/5h00h1262182684.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 = 61
Frequency = 1
1 2 3 4 5 6
-1306.71723 -1755.87107 -1212.73867 -1350.21064 -1002.68901 435.31859
7 8 9 10 11 12
-1183.03853 -904.08753 47.34052 166.62171 -588.49885 -682.23179
13 14 15 16 17 18
-98.23367 -359.28818 551.21561 -383.13967 1073.80404 1388.83722
19 20 21 22 23 24
1008.97978 739.62659 1297.76619 700.85432 1087.61136 312.34969
25 26 27 28 29 30
643.67413 401.61371 590.62082 -413.10779 562.73468 -307.26508
31 32 33 34 35 36
417.20032 -142.68813 258.09901 -829.71120 35.68760 166.87945
37 38 39 40 41 42
852.59671 689.84647 53.91693 935.31679 715.84598 1321.85752
43 44 45 46 47 48
1224.99594 2675.43777 990.38408 1474.00714 174.06287 928.05616
49 50 51 52 53 54
490.94072 1023.69907 16.98531 1211.14131 -1349.69570 -2838.74825
55 56 57 58 59 60
-1468.13751 -2368.28870 -2593.58981 -1511.77197 -708.86299 -725.05350
61
-582.26067
> postscript(file="/var/www/html/rcomp/tmp/6uohg1262182684.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 = 61
Frequency = 1
lag(myerror, k = 1) myerror
0 -1306.71723 NA
1 -1755.87107 -1306.71723
2 -1212.73867 -1755.87107
3 -1350.21064 -1212.73867
4 -1002.68901 -1350.21064
5 435.31859 -1002.68901
6 -1183.03853 435.31859
7 -904.08753 -1183.03853
8 47.34052 -904.08753
9 166.62171 47.34052
10 -588.49885 166.62171
11 -682.23179 -588.49885
12 -98.23367 -682.23179
13 -359.28818 -98.23367
14 551.21561 -359.28818
15 -383.13967 551.21561
16 1073.80404 -383.13967
17 1388.83722 1073.80404
18 1008.97978 1388.83722
19 739.62659 1008.97978
20 1297.76619 739.62659
21 700.85432 1297.76619
22 1087.61136 700.85432
23 312.34969 1087.61136
24 643.67413 312.34969
25 401.61371 643.67413
26 590.62082 401.61371
27 -413.10779 590.62082
28 562.73468 -413.10779
29 -307.26508 562.73468
30 417.20032 -307.26508
31 -142.68813 417.20032
32 258.09901 -142.68813
33 -829.71120 258.09901
34 35.68760 -829.71120
35 166.87945 35.68760
36 852.59671 166.87945
37 689.84647 852.59671
38 53.91693 689.84647
39 935.31679 53.91693
40 715.84598 935.31679
41 1321.85752 715.84598
42 1224.99594 1321.85752
43 2675.43777 1224.99594
44 990.38408 2675.43777
45 1474.00714 990.38408
46 174.06287 1474.00714
47 928.05616 174.06287
48 490.94072 928.05616
49 1023.69907 490.94072
50 16.98531 1023.69907
51 1211.14131 16.98531
52 -1349.69570 1211.14131
53 -2838.74825 -1349.69570
54 -1468.13751 -2838.74825
55 -2368.28870 -1468.13751
56 -2593.58981 -2368.28870
57 -1511.77197 -2593.58981
58 -708.86299 -1511.77197
59 -725.05350 -708.86299
60 -582.26067 -725.05350
61 NA -582.26067
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -1755.87107 -1306.71723
[2,] -1212.73867 -1755.87107
[3,] -1350.21064 -1212.73867
[4,] -1002.68901 -1350.21064
[5,] 435.31859 -1002.68901
[6,] -1183.03853 435.31859
[7,] -904.08753 -1183.03853
[8,] 47.34052 -904.08753
[9,] 166.62171 47.34052
[10,] -588.49885 166.62171
[11,] -682.23179 -588.49885
[12,] -98.23367 -682.23179
[13,] -359.28818 -98.23367
[14,] 551.21561 -359.28818
[15,] -383.13967 551.21561
[16,] 1073.80404 -383.13967
[17,] 1388.83722 1073.80404
[18,] 1008.97978 1388.83722
[19,] 739.62659 1008.97978
[20,] 1297.76619 739.62659
[21,] 700.85432 1297.76619
[22,] 1087.61136 700.85432
[23,] 312.34969 1087.61136
[24,] 643.67413 312.34969
[25,] 401.61371 643.67413
[26,] 590.62082 401.61371
[27,] -413.10779 590.62082
[28,] 562.73468 -413.10779
[29,] -307.26508 562.73468
[30,] 417.20032 -307.26508
[31,] -142.68813 417.20032
[32,] 258.09901 -142.68813
[33,] -829.71120 258.09901
[34,] 35.68760 -829.71120
[35,] 166.87945 35.68760
[36,] 852.59671 166.87945
[37,] 689.84647 852.59671
[38,] 53.91693 689.84647
[39,] 935.31679 53.91693
[40,] 715.84598 935.31679
[41,] 1321.85752 715.84598
[42,] 1224.99594 1321.85752
[43,] 2675.43777 1224.99594
[44,] 990.38408 2675.43777
[45,] 1474.00714 990.38408
[46,] 174.06287 1474.00714
[47,] 928.05616 174.06287
[48,] 490.94072 928.05616
[49,] 1023.69907 490.94072
[50,] 16.98531 1023.69907
[51,] 1211.14131 16.98531
[52,] -1349.69570 1211.14131
[53,] -2838.74825 -1349.69570
[54,] -1468.13751 -2838.74825
[55,] -2368.28870 -1468.13751
[56,] -2593.58981 -2368.28870
[57,] -1511.77197 -2593.58981
[58,] -708.86299 -1511.77197
[59,] -725.05350 -708.86299
[60,] -582.26067 -725.05350
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -1755.87107 -1306.71723
2 -1212.73867 -1755.87107
3 -1350.21064 -1212.73867
4 -1002.68901 -1350.21064
5 435.31859 -1002.68901
6 -1183.03853 435.31859
7 -904.08753 -1183.03853
8 47.34052 -904.08753
9 166.62171 47.34052
10 -588.49885 166.62171
11 -682.23179 -588.49885
12 -98.23367 -682.23179
13 -359.28818 -98.23367
14 551.21561 -359.28818
15 -383.13967 551.21561
16 1073.80404 -383.13967
17 1388.83722 1073.80404
18 1008.97978 1388.83722
19 739.62659 1008.97978
20 1297.76619 739.62659
21 700.85432 1297.76619
22 1087.61136 700.85432
23 312.34969 1087.61136
24 643.67413 312.34969
25 401.61371 643.67413
26 590.62082 401.61371
27 -413.10779 590.62082
28 562.73468 -413.10779
29 -307.26508 562.73468
30 417.20032 -307.26508
31 -142.68813 417.20032
32 258.09901 -142.68813
33 -829.71120 258.09901
34 35.68760 -829.71120
35 166.87945 35.68760
36 852.59671 166.87945
37 689.84647 852.59671
38 53.91693 689.84647
39 935.31679 53.91693
40 715.84598 935.31679
41 1321.85752 715.84598
42 1224.99594 1321.85752
43 2675.43777 1224.99594
44 990.38408 2675.43777
45 1474.00714 990.38408
46 174.06287 1474.00714
47 928.05616 174.06287
48 490.94072 928.05616
49 1023.69907 490.94072
50 16.98531 1023.69907
51 1211.14131 16.98531
52 -1349.69570 1211.14131
53 -2838.74825 -1349.69570
54 -1468.13751 -2838.74825
55 -2368.28870 -1468.13751
56 -2593.58981 -2368.28870
57 -1511.77197 -2593.58981
58 -708.86299 -1511.77197
59 -725.05350 -708.86299
60 -582.26067 -725.05350
> 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/7xdsk1262182684.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/84a011262182684.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/9vk0s1262182684.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/10vynp1262182684.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/11ke5d1262182684.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/125vkz1262182684.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/13ydwm1262182685.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/14pmqc1262182685.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/15chce1262182685.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/16nn9w1262182685.tab")
+ }
>
> try(system("convert tmp/19cz91262182684.ps tmp/19cz91262182684.png",intern=TRUE))
character(0)
> try(system("convert tmp/26d5x1262182684.ps tmp/26d5x1262182684.png",intern=TRUE))
character(0)
> try(system("convert tmp/30i061262182684.ps tmp/30i061262182684.png",intern=TRUE))
character(0)
> try(system("convert tmp/4qb2a1262182684.ps tmp/4qb2a1262182684.png",intern=TRUE))
character(0)
> try(system("convert tmp/5h00h1262182684.ps tmp/5h00h1262182684.png",intern=TRUE))
character(0)
> try(system("convert tmp/6uohg1262182684.ps tmp/6uohg1262182684.png",intern=TRUE))
character(0)
> try(system("convert tmp/7xdsk1262182684.ps tmp/7xdsk1262182684.png",intern=TRUE))
character(0)
> try(system("convert tmp/84a011262182684.ps tmp/84a011262182684.png",intern=TRUE))
character(0)
> try(system("convert tmp/9vk0s1262182684.ps tmp/9vk0s1262182684.png",intern=TRUE))
character(0)
> try(system("convert tmp/10vynp1262182684.ps tmp/10vynp1262182684.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.374 1.569 3.465