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 = 'No 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
1 10519.2 1154.8 1 0 0 0 0 0 0 0 0 0 0
2 10414.9 1206.7 0 1 0 0 0 0 0 0 0 0 0
3 12476.8 1199.0 0 0 1 0 0 0 0 0 0 0 0
4 12384.6 1265.0 0 0 0 1 0 0 0 0 0 0 0
5 12266.7 1247.1 0 0 0 0 1 0 0 0 0 0 0
6 12919.9 1116.5 0 0 0 0 0 1 0 0 0 0 0
7 11497.3 1153.9 0 0 0 0 0 0 1 0 0 0 0
8 12142.0 1077.4 0 0 0 0 0 0 0 1 0 0 0
9 13919.4 1132.5 0 0 0 0 0 0 0 0 1 0 0
10 12656.8 1058.8 0 0 0 0 0 0 0 0 0 1 0
11 12034.1 1195.1 0 0 0 0 0 0 0 0 0 0 1
12 13199.7 1263.4 0 0 0 0 0 0 0 0 0 0 0
13 10881.3 1023.1 1 0 0 0 0 0 0 0 0 0 0
14 11301.2 1141.0 0 1 0 0 0 0 0 0 0 0 0
15 13643.9 1116.3 0 0 1 0 0 0 0 0 0 0 0
16 12517.0 1135.6 0 0 0 1 0 0 0 0 0 0 0
17 13981.1 1210.5 0 0 0 0 1 0 0 0 0 0 0
18 14275.7 1230.0 0 0 0 0 0 1 0 0 0 0 0
19 13425.0 1136.5 0 0 0 0 0 0 1 0 0 0 0
20 13565.7 1068.7 0 0 0 0 0 0 0 1 0 0 0
21 16216.3 1372.5 0 0 0 0 0 0 0 0 1 0 0
22 12970.0 1049.9 0 0 0 0 0 0 0 0 0 1 0
23 14079.9 1302.2 0 0 0 0 0 0 0 0 0 0 1
24 14235.0 1305.9 0 0 0 0 0 0 0 0 0 0 0
25 12213.4 1173.5 1 0 0 0 0 0 0 0 0 0 0
26 12581.0 1277.4 0 1 0 0 0 0 0 0 0 0 0
27 14130.4 1238.6 0 0 1 0 0 0 0 0 0 0 0
28 14210.8 1508.6 0 0 0 1 0 0 0 0 0 0 0
29 14378.5 1423.4 0 0 0 0 1 0 0 0 0 0 0
30 13142.8 1375.1 0 0 0 0 0 1 0 0 0 0 0
31 13714.7 1344.1 0 0 0 0 0 0 1 0 0 0 0
32 13621.9 1287.5 0 0 0 0 0 0 0 1 0 0 0
33 15379.8 1446.9 0 0 0 0 0 0 0 0 1 0 0
34 13306.3 1451.0 0 0 0 0 0 0 0 0 0 1 0
35 14391.2 1604.4 0 0 0 0 0 0 0 0 0 0 1
36 14909.9 1501.5 0 0 0 0 0 0 0 0 0 0 0
37 14025.4 1522.8 1 0 0 0 0 0 0 0 0 0 0
38 12951.2 1328.0 0 1 0 0 0 0 0 0 0 0 0
39 14344.3 1420.5 0 0 1 0 0 0 0 0 0 0 0
40 16093.4 1648.0 0 0 0 1 0 0 0 0 0 0 0
41 15413.6 1631.1 0 0 0 0 1 0 0 0 0 0 0
42 14705.7 1396.6 0 0 0 0 0 1 0 0 0 0 0
43 15972.8 1663.4 0 0 0 0 0 0 1 0 0 0 0
44 16241.4 1283.0 0 0 0 0 0 0 0 1 0 0 0
45 16626.4 1582.4 0 0 0 0 0 0 0 0 1 0 0
46 17136.2 1785.2 0 0 0 0 0 0 0 0 0 1 0
47 15622.9 1853.6 0 0 0 0 0 0 0 0 0 0 1
48 18003.9 1994.1 0 0 0 0 0 0 0 0 0 0 0
49 16136.1 2042.8 1 0 0 0 0 0 0 0 0 0 0
50 14423.7 1586.1 0 1 0 0 0 0 0 0 0 0 0
51 16789.4 1942.4 0 0 1 0 0 0 0 0 0 0 0
52 16782.2 1763.6 0 0 0 1 0 0 0 0 0 0 0
53 14133.8 1819.9 0 0 0 0 1 0 0 0 0 0 0
54 12607.0 1836.0 0 0 0 0 0 1 0 0 0 0 0
55 12004.5 1447.5 0 0 0 0 0 0 1 0 0 0 0
56 12175.4 1509.5 0 0 0 0 0 0 0 1 0 0 0
57 13268.0 1661.2 0 0 0 0 0 0 0 0 1 0 0
58 12299.3 1456.2 0 0 0 0 0 0 0 0 0 1 0
59 11800.6 1310.9 0 0 0 0 0 0 0 0 0 0 1
60 13873.3 1542.1 0 0 0 0 0 0 0 0 0 0 0
61 12315.0 1537.7 1 0 0 0 0 0 0 0 0 0 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) InvoerAM M1 M2 M3 M4
8364.552 4.259 -1684.400 -1600.385 20.527 -202.969
M5 M6 M7 M8 M9 M10
-575.369 -758.071 -787.570 -118.799 588.145 -484.157
M11
-968.319
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-2819.2 -754.1 264.6 768.2 2531.2
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8364.5521 1129.8583 7.403 1.78e-09 ***
InvoerAM 4.2591 0.6485 6.568 3.38e-08 ***
M1 -1684.3997 749.1309 -2.248 0.0292 *
M2 -1600.3847 790.9554 -2.023 0.0486 *
M3 20.5274 783.8650 0.026 0.9792
M4 -202.9686 779.6212 -0.260 0.7957
M5 -575.3690 779.5533 -0.738 0.4641
M6 -758.0708 783.3259 -0.968 0.3380
M7 -787.5704 786.7137 -1.001 0.3218
M8 -118.7994 799.0662 -0.149 0.8824
M9 588.1446 780.5637 0.753 0.4548
M10 -484.1569 785.7201 -0.616 0.5407
M11 -968.3192 779.9904 -1.241 0.2205
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1231 on 48 degrees of freedom
Multiple R-squared: 0.5849, Adjusted R-squared: 0.4811
F-statistic: 5.636 on 12 and 48 DF, p-value: 6.459e-06
> 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,] 4.863360e-02 9.726720e-02 0.9513664
[2,] 8.549180e-02 1.709836e-01 0.9145082
[3,] 1.363204e-01 2.726407e-01 0.8636796
[4,] 1.682960e-01 3.365919e-01 0.8317040
[5,] 1.420051e-01 2.840102e-01 0.8579949
[6,] 1.940377e-01 3.880754e-01 0.8059623
[7,] 1.256603e-01 2.513205e-01 0.8743397
[8,] 1.215754e-01 2.431509e-01 0.8784246
[9,] 8.183646e-02 1.636729e-01 0.9181635
[10,] 6.109964e-02 1.221993e-01 0.9389004
[11,] 4.506834e-02 9.013668e-02 0.9549317
[12,] 2.690048e-02 5.380095e-02 0.9730995
[13,] 1.647033e-02 3.294065e-02 0.9835297
[14,] 8.897277e-03 1.779455e-02 0.9911027
[15,] 9.231151e-03 1.846230e-02 0.9907688
[16,] 4.767255e-03 9.534510e-03 0.9952327
[17,] 2.324075e-03 4.648151e-03 0.9976759
[18,] 1.364506e-03 2.729012e-03 0.9986355
[19,] 9.586987e-04 1.917397e-03 0.9990413
[20,] 4.155872e-04 8.311745e-04 0.9995844
[21,] 1.716429e-04 3.432857e-04 0.9998284
[22,] 1.146950e-04 2.293899e-04 0.9998853
[23,] 5.857986e-05 1.171597e-04 0.9999414
[24,] 2.025303e-05 4.050607e-05 0.9999797
[25,] 1.307917e-05 2.615835e-05 0.9999869
[26,] 6.774786e-06 1.354957e-05 0.9999932
[27,] 3.983890e-05 7.967780e-05 0.9999602
[28,] 3.311950e-05 6.623900e-05 0.9999669
[29,] 1.595389e-02 3.190777e-02 0.9840461
[30,] 3.321526e-01 6.643051e-01 0.6678474
> postscript(file="/var/www/html/rcomp/tmp/1pbb31262181729.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/2sihu1262181729.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/30na61262181729.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/4lvvf1262181729.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/510jm1262181729.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 7
-1079.3710 -1488.7337 -1014.9507 -1164.7559 -834.0175 558.1239 -994.2672
8 9 10 11 12 13 14
-692.5164 143.2627 266.8605 -452.1937 -545.8100 -156.3464 -322.6103
15 16 17 18 19 20 21
504.3776 -481.2273 1036.2659 1430.5150 1007.5413 768.2379 1417.9766
22 23 24 25 26 27 28
617.9666 1137.4558 308.4779 535.1837 376.2473 469.9886 -376.0748
29 30 31 32 33 34 35
526.9017 -320.3816 413.0504 -107.4551 264.5990 -754.0619 161.6531
36 37 38 39 40 41 42
150.2963 859.4770 530.9364 -90.8433 912.8055 677.3848 1150.9475
43 44 45 46 47 48 49
1311.2170 2531.2109 934.0897 1652.4440 331.9832 1146.2593 755.4405
50 51 52 53 54 55 56
904.1604 131.4279 1109.2525 -1406.5349 -2819.2048 -1737.5415 -2499.4772
57 58 59 60 61
-2759.9280 -1783.2093 -1178.8985 -1059.2235 -914.3837
> postscript(file="/var/www/html/rcomp/tmp/6jvi81262181729.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 -1079.3710 NA
1 -1488.7337 -1079.3710
2 -1014.9507 -1488.7337
3 -1164.7559 -1014.9507
4 -834.0175 -1164.7559
5 558.1239 -834.0175
6 -994.2672 558.1239
7 -692.5164 -994.2672
8 143.2627 -692.5164
9 266.8605 143.2627
10 -452.1937 266.8605
11 -545.8100 -452.1937
12 -156.3464 -545.8100
13 -322.6103 -156.3464
14 504.3776 -322.6103
15 -481.2273 504.3776
16 1036.2659 -481.2273
17 1430.5150 1036.2659
18 1007.5413 1430.5150
19 768.2379 1007.5413
20 1417.9766 768.2379
21 617.9666 1417.9766
22 1137.4558 617.9666
23 308.4779 1137.4558
24 535.1837 308.4779
25 376.2473 535.1837
26 469.9886 376.2473
27 -376.0748 469.9886
28 526.9017 -376.0748
29 -320.3816 526.9017
30 413.0504 -320.3816
31 -107.4551 413.0504
32 264.5990 -107.4551
33 -754.0619 264.5990
34 161.6531 -754.0619
35 150.2963 161.6531
36 859.4770 150.2963
37 530.9364 859.4770
38 -90.8433 530.9364
39 912.8055 -90.8433
40 677.3848 912.8055
41 1150.9475 677.3848
42 1311.2170 1150.9475
43 2531.2109 1311.2170
44 934.0897 2531.2109
45 1652.4440 934.0897
46 331.9832 1652.4440
47 1146.2593 331.9832
48 755.4405 1146.2593
49 904.1604 755.4405
50 131.4279 904.1604
51 1109.2525 131.4279
52 -1406.5349 1109.2525
53 -2819.2048 -1406.5349
54 -1737.5415 -2819.2048
55 -2499.4772 -1737.5415
56 -2759.9280 -2499.4772
57 -1783.2093 -2759.9280
58 -1178.8985 -1783.2093
59 -1059.2235 -1178.8985
60 -914.3837 -1059.2235
61 NA -914.3837
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -1488.7337 -1079.3710
[2,] -1014.9507 -1488.7337
[3,] -1164.7559 -1014.9507
[4,] -834.0175 -1164.7559
[5,] 558.1239 -834.0175
[6,] -994.2672 558.1239
[7,] -692.5164 -994.2672
[8,] 143.2627 -692.5164
[9,] 266.8605 143.2627
[10,] -452.1937 266.8605
[11,] -545.8100 -452.1937
[12,] -156.3464 -545.8100
[13,] -322.6103 -156.3464
[14,] 504.3776 -322.6103
[15,] -481.2273 504.3776
[16,] 1036.2659 -481.2273
[17,] 1430.5150 1036.2659
[18,] 1007.5413 1430.5150
[19,] 768.2379 1007.5413
[20,] 1417.9766 768.2379
[21,] 617.9666 1417.9766
[22,] 1137.4558 617.9666
[23,] 308.4779 1137.4558
[24,] 535.1837 308.4779
[25,] 376.2473 535.1837
[26,] 469.9886 376.2473
[27,] -376.0748 469.9886
[28,] 526.9017 -376.0748
[29,] -320.3816 526.9017
[30,] 413.0504 -320.3816
[31,] -107.4551 413.0504
[32,] 264.5990 -107.4551
[33,] -754.0619 264.5990
[34,] 161.6531 -754.0619
[35,] 150.2963 161.6531
[36,] 859.4770 150.2963
[37,] 530.9364 859.4770
[38,] -90.8433 530.9364
[39,] 912.8055 -90.8433
[40,] 677.3848 912.8055
[41,] 1150.9475 677.3848
[42,] 1311.2170 1150.9475
[43,] 2531.2109 1311.2170
[44,] 934.0897 2531.2109
[45,] 1652.4440 934.0897
[46,] 331.9832 1652.4440
[47,] 1146.2593 331.9832
[48,] 755.4405 1146.2593
[49,] 904.1604 755.4405
[50,] 131.4279 904.1604
[51,] 1109.2525 131.4279
[52,] -1406.5349 1109.2525
[53,] -2819.2048 -1406.5349
[54,] -1737.5415 -2819.2048
[55,] -2499.4772 -1737.5415
[56,] -2759.9280 -2499.4772
[57,] -1783.2093 -2759.9280
[58,] -1178.8985 -1783.2093
[59,] -1059.2235 -1178.8985
[60,] -914.3837 -1059.2235
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -1488.7337 -1079.3710
2 -1014.9507 -1488.7337
3 -1164.7559 -1014.9507
4 -834.0175 -1164.7559
5 558.1239 -834.0175
6 -994.2672 558.1239
7 -692.5164 -994.2672
8 143.2627 -692.5164
9 266.8605 143.2627
10 -452.1937 266.8605
11 -545.8100 -452.1937
12 -156.3464 -545.8100
13 -322.6103 -156.3464
14 504.3776 -322.6103
15 -481.2273 504.3776
16 1036.2659 -481.2273
17 1430.5150 1036.2659
18 1007.5413 1430.5150
19 768.2379 1007.5413
20 1417.9766 768.2379
21 617.9666 1417.9766
22 1137.4558 617.9666
23 308.4779 1137.4558
24 535.1837 308.4779
25 376.2473 535.1837
26 469.9886 376.2473
27 -376.0748 469.9886
28 526.9017 -376.0748
29 -320.3816 526.9017
30 413.0504 -320.3816
31 -107.4551 413.0504
32 264.5990 -107.4551
33 -754.0619 264.5990
34 161.6531 -754.0619
35 150.2963 161.6531
36 859.4770 150.2963
37 530.9364 859.4770
38 -90.8433 530.9364
39 912.8055 -90.8433
40 677.3848 912.8055
41 1150.9475 677.3848
42 1311.2170 1150.9475
43 2531.2109 1311.2170
44 934.0897 2531.2109
45 1652.4440 934.0897
46 331.9832 1652.4440
47 1146.2593 331.9832
48 755.4405 1146.2593
49 904.1604 755.4405
50 131.4279 904.1604
51 1109.2525 131.4279
52 -1406.5349 1109.2525
53 -2819.2048 -1406.5349
54 -1737.5415 -2819.2048
55 -2499.4772 -1737.5415
56 -2759.9280 -2499.4772
57 -1783.2093 -2759.9280
58 -1178.8985 -1783.2093
59 -1059.2235 -1178.8985
60 -914.3837 -1059.2235
> 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/7qgka1262181729.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/8fd8z1262181729.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/97y7a1262181729.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/10l8pk1262181729.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/11bx9f1262181729.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/126t2h1262181730.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/13n64e1262181730.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/14e5t51262181730.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/15ywkg1262181730.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/16xcwn1262181730.tab")
+ }
>
> try(system("convert tmp/1pbb31262181729.ps tmp/1pbb31262181729.png",intern=TRUE))
character(0)
> try(system("convert tmp/2sihu1262181729.ps tmp/2sihu1262181729.png",intern=TRUE))
character(0)
> try(system("convert tmp/30na61262181729.ps tmp/30na61262181729.png",intern=TRUE))
character(0)
> try(system("convert tmp/4lvvf1262181729.ps tmp/4lvvf1262181729.png",intern=TRUE))
character(0)
> try(system("convert tmp/510jm1262181729.ps tmp/510jm1262181729.png",intern=TRUE))
character(0)
> try(system("convert tmp/6jvi81262181729.ps tmp/6jvi81262181729.png",intern=TRUE))
character(0)
> try(system("convert tmp/7qgka1262181729.ps tmp/7qgka1262181729.png",intern=TRUE))
character(0)
> try(system("convert tmp/8fd8z1262181729.ps tmp/8fd8z1262181729.png",intern=TRUE))
character(0)
> try(system("convert tmp/97y7a1262181729.ps tmp/97y7a1262181729.png",intern=TRUE))
character(0)
> try(system("convert tmp/10l8pk1262181729.ps tmp/10l8pk1262181729.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.440 1.556 3.788