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(427.25
+ ,1113.89
+ ,1144.94
+ ,1131.13
+ ,1111.92
+ ,391.25
+ ,1107.3
+ ,1113.89
+ ,1144.94
+ ,1131.13
+ ,397.20
+ ,1120.68
+ ,1107.3
+ ,1113.89
+ ,1144.94
+ ,394.80
+ ,1140.84
+ ,1120.68
+ ,1107.3
+ ,1113.89
+ ,391.50
+ ,1101.72
+ ,1140.84
+ ,1120.68
+ ,1107.3
+ ,407.65
+ ,1104.24
+ ,1101.72
+ ,1140.84
+ ,1120.68
+ ,418.10
+ ,1114.58
+ ,1104.24
+ ,1101.72
+ ,1140.84
+ ,429.10
+ ,1130.2
+ ,1114.58
+ ,1104.24
+ ,1101.72
+ ,452.85
+ ,1173.78
+ ,1130.2
+ ,1114.58
+ ,1104.24
+ ,427.75
+ ,1211.92
+ ,1173.78
+ ,1130.2
+ ,1114.58
+ ,420.90
+ ,1181.27
+ ,1211.92
+ ,1173.78
+ ,1130.2
+ ,433.45
+ ,1203.6
+ ,1181.27
+ ,1211.92
+ ,1173.78
+ ,427.15
+ ,1180.59
+ ,1203.6
+ ,1181.27
+ ,1211.92
+ ,427.90
+ ,1156.85
+ ,1180.59
+ ,1203.6
+ ,1181.27
+ ,415.35
+ ,1191.5
+ ,1156.85
+ ,1180.59
+ ,1203.6
+ ,432.60
+ ,1191.33
+ ,1191.5
+ ,1156.85
+ ,1180.59
+ ,431.65
+ ,1234.18
+ ,1191.33
+ ,1191.5
+ ,1156.85
+ ,439.60
+ ,1220.33
+ ,1234.18
+ ,1191.33
+ ,1191.5
+ ,466.10
+ ,1228.81
+ ,1220.33
+ ,1234.18
+ ,1191.33
+ ,459.50
+ ,1207.01
+ ,1228.81
+ ,1220.33
+ ,1234.18
+ ,499.75
+ ,1249.48
+ ,1207.01
+ ,1228.81
+ ,1220.33
+ ,530.00
+ ,1248.29
+ ,1249.48
+ ,1207.01
+ ,1228.81
+ ,568.25
+ ,1280.08
+ ,1248.29
+ ,1249.48
+ ,1207.01
+ ,564.25
+ ,1280.66
+ ,1280.08
+ ,1248.29
+ ,1249.48
+ ,587.00
+ ,1302.88
+ ,1280.66
+ ,1280.08
+ ,1248.29
+ ,661.00
+ ,1310.61
+ ,1302.88
+ ,1280.66
+ ,1280.08
+ ,625.00
+ ,1270.05
+ ,1310.61
+ ,1302.88
+ ,1280.66
+ ,622.95
+ ,1270.06
+ ,1270.05
+ ,1310.61
+ ,1302.88
+ ,637.25
+ ,1278.53
+ ,1270.06
+ ,1270.05
+ ,1310.61
+ ,621.05
+ ,1303.8
+ ,1278.53
+ ,1270.06
+ ,1270.05
+ ,600.60
+ ,1335.83
+ ,1303.8
+ ,1278.53
+ ,1270.06
+ ,614.10
+ ,1377.76
+ ,1335.83
+ ,1303.8
+ ,1278.53
+ ,648.75
+ ,1400.63
+ ,1377.76
+ ,1335.83
+ ,1303.8
+ ,639.75
+ ,1418.03
+ ,1400.63
+ ,1377.76
+ ,1335.83
+ ,660.20
+ ,1437.9
+ ,1418.03
+ ,1400.63
+ ,1377.76
+ ,670.40
+ ,1406.8
+ ,1437.9
+ ,1418.03
+ ,1400.63
+ ,658.25
+ ,1420.83
+ ,1406.8
+ ,1437.9
+ ,1418.03
+ ,673.60
+ ,1482.37
+ ,1420.83
+ ,1406.8
+ ,1437.9
+ ,666.50
+ ,1530.63
+ ,1482.37
+ ,1420.83
+ ,1406.8
+ ,654.75
+ ,1504.66
+ ,1530.63
+ ,1482.37
+ ,1420.83
+ ,665.75
+ ,1455.18
+ ,1504.66
+ ,1530.63
+ ,1482.37
+ ,672.00
+ ,1473.96
+ ,1455.18
+ ,1504.66
+ ,1530.63
+ ,742.50
+ ,1527.29
+ ,1473.96
+ ,1455.18
+ ,1504.66
+ ,790.25
+ ,1545.79
+ ,1527.29
+ ,1473.96
+ ,1455.18
+ ,784.25
+ ,1479.63
+ ,1545.79
+ ,1527.29
+ ,1473.96
+ ,846.75
+ ,1467.97
+ ,1479.63
+ ,1545.79
+ ,1527.29
+ ,914.75
+ ,1378.6
+ ,1467.97
+ ,1479.63
+ ,1545.79
+ ,988.50
+ ,1330.45
+ ,1378.6
+ ,1467.97
+ ,1479.63
+ ,887.75
+ ,1326.41
+ ,1330.45
+ ,1378.6
+ ,1467.97
+ ,853.00
+ ,1385.97
+ ,1326.41
+ ,1330.45
+ ,1378.6
+ ,888.25
+ ,1399.62
+ ,1385.97
+ ,1326.41
+ ,1330.45
+ ,937.50
+ ,1276.69
+ ,1399.62
+ ,1385.97
+ ,1326.41
+ ,912.50
+ ,1269.42
+ ,1276.69
+ ,1399.62
+ ,1385.97
+ ,822.25
+ ,1287.83
+ ,1269.42
+ ,1276.69
+ ,1399.62
+ ,880.00
+ ,1164.17
+ ,1287.83
+ ,1269.42
+ ,1276.69
+ ,729.50
+ ,968.67
+ ,1164.17
+ ,1287.83
+ ,1269.42
+ ,778.00
+ ,888.61
+ ,968.67
+ ,1164.17
+ ,1287.83)
+ ,dim=c(5
+ ,57)
+ ,dimnames=list(c('x(t)'
+ ,'y(t)'
+ ,'y(t-1)'
+ ,'y(t-2)'
+ ,'y(t-3)
')
+ ,1:57))
> y <- array(NA,dim=c(5,57),dimnames=list(c('x(t)','y(t)','y(t-1)','y(t-2)','y(t-3)
'),1:57))
> 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 = '2'
> #'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
y(t) x(t) y(t-1) y(t-2) y(t-3)\r
1 1113.89 427.25 1144.94 1131.13 1111.92
2 1107.30 391.25 1113.89 1144.94 1131.13
3 1120.68 397.20 1107.30 1113.89 1144.94
4 1140.84 394.80 1120.68 1107.30 1113.89
5 1101.72 391.50 1140.84 1120.68 1107.30
6 1104.24 407.65 1101.72 1140.84 1120.68
7 1114.58 418.10 1104.24 1101.72 1140.84
8 1130.20 429.10 1114.58 1104.24 1101.72
9 1173.78 452.85 1130.20 1114.58 1104.24
10 1211.92 427.75 1173.78 1130.20 1114.58
11 1181.27 420.90 1211.92 1173.78 1130.20
12 1203.60 433.45 1181.27 1211.92 1173.78
13 1180.59 427.15 1203.60 1181.27 1211.92
14 1156.85 427.90 1180.59 1203.60 1181.27
15 1191.50 415.35 1156.85 1180.59 1203.60
16 1191.33 432.60 1191.50 1156.85 1180.59
17 1234.18 431.65 1191.33 1191.50 1156.85
18 1220.33 439.60 1234.18 1191.33 1191.50
19 1228.81 466.10 1220.33 1234.18 1191.33
20 1207.01 459.50 1228.81 1220.33 1234.18
21 1249.48 499.75 1207.01 1228.81 1220.33
22 1248.29 530.00 1249.48 1207.01 1228.81
23 1280.08 568.25 1248.29 1249.48 1207.01
24 1280.66 564.25 1280.08 1248.29 1249.48
25 1302.88 587.00 1280.66 1280.08 1248.29
26 1310.61 661.00 1302.88 1280.66 1280.08
27 1270.05 625.00 1310.61 1302.88 1280.66
28 1270.06 622.95 1270.05 1310.61 1302.88
29 1278.53 637.25 1270.06 1270.05 1310.61
30 1303.80 621.05 1278.53 1270.06 1270.05
31 1335.83 600.60 1303.80 1278.53 1270.06
32 1377.76 614.10 1335.83 1303.80 1278.53
33 1400.63 648.75 1377.76 1335.83 1303.80
34 1418.03 639.75 1400.63 1377.76 1335.83
35 1437.90 660.20 1418.03 1400.63 1377.76
36 1406.80 670.40 1437.90 1418.03 1400.63
37 1420.83 658.25 1406.80 1437.90 1418.03
38 1482.37 673.60 1420.83 1406.80 1437.90
39 1530.63 666.50 1482.37 1420.83 1406.80
40 1504.66 654.75 1530.63 1482.37 1420.83
41 1455.18 665.75 1504.66 1530.63 1482.37
42 1473.96 672.00 1455.18 1504.66 1530.63
43 1527.29 742.50 1473.96 1455.18 1504.66
44 1545.79 790.25 1527.29 1473.96 1455.18
45 1479.63 784.25 1545.79 1527.29 1473.96
46 1467.97 846.75 1479.63 1545.79 1527.29
47 1378.60 914.75 1467.97 1479.63 1545.79
48 1330.45 988.50 1378.60 1467.97 1479.63
49 1326.41 887.75 1330.45 1378.60 1467.97
50 1385.97 853.00 1326.41 1330.45 1378.60
51 1399.62 888.25 1385.97 1326.41 1330.45
52 1276.69 937.50 1399.62 1385.97 1326.41
53 1269.42 912.50 1276.69 1399.62 1385.97
54 1287.83 822.25 1269.42 1276.69 1399.62
55 1164.17 880.00 1287.83 1269.42 1276.69
56 968.67 729.50 1164.17 1287.83 1269.42
57 888.61 778.00 968.67 1164.17 1287.83
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) `x(t)` `y(t-1)` `y(t-2)` `y(t-3)\r`
-101.3731 -0.1770 1.3382 -0.6129 0.4352
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-121.87 -30.63 13.28 24.84 78.78
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -101.3731 72.0886 -1.406 0.16560
`x(t)` -0.1770 0.0585 -3.026 0.00384 **
`y(t-1)` 1.3382 0.1290 10.370 2.92e-14 ***
`y(t-2)` -0.6129 0.2234 -2.743 0.00832 **
`y(t-3)\r` 0.4352 0.1684 2.584 0.01262 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 41.94 on 52 degrees of freedom
Multiple R-squared: 0.9201, Adjusted R-squared: 0.914
F-statistic: 149.7 on 4 and 52 DF, p-value: < 2.2e-16
> 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.459865e-02 8.919731e-02 0.9554013
[2,] 8.458051e-02 1.691610e-01 0.9154195
[3,] 1.584147e-01 3.168295e-01 0.8415853
[4,] 9.671541e-02 1.934308e-01 0.9032846
[5,] 6.087304e-02 1.217461e-01 0.9391270
[6,] 5.436680e-02 1.087336e-01 0.9456332
[7,] 3.468273e-02 6.936546e-02 0.9653173
[8,] 3.345703e-02 6.691405e-02 0.9665430
[9,] 1.876289e-02 3.752578e-02 0.9812371
[10,] 2.745199e-02 5.490399e-02 0.9725480
[11,] 1.777094e-02 3.554189e-02 0.9822291
[12,] 9.997847e-03 1.999569e-02 0.9900022
[13,] 1.038609e-02 2.077218e-02 0.9896139
[14,] 5.821952e-03 1.164390e-02 0.9941780
[15,] 5.473164e-03 1.094633e-02 0.9945268
[16,] 3.571044e-03 7.142088e-03 0.9964290
[17,] 2.234110e-03 4.468220e-03 0.9977659
[18,] 1.227715e-03 2.455430e-03 0.9987723
[19,] 8.668059e-04 1.733612e-03 0.9991332
[20,] 2.049361e-03 4.098723e-03 0.9979506
[21,] 1.211025e-03 2.422050e-03 0.9987890
[22,] 5.945165e-04 1.189033e-03 0.9994055
[23,] 3.296831e-04 6.593662e-04 0.9996703
[24,] 2.581280e-04 5.162559e-04 0.9997419
[25,] 3.337164e-04 6.674328e-04 0.9996663
[26,] 2.080420e-04 4.160839e-04 0.9997920
[27,] 1.316149e-04 2.632299e-04 0.9998684
[28,] 7.720245e-05 1.544049e-04 0.9999228
[29,] 6.029629e-05 1.205926e-04 0.9999397
[30,] 3.609290e-05 7.218580e-05 0.9999639
[31,] 9.201596e-05 1.840319e-04 0.9999080
[32,] 1.252606e-04 2.505211e-04 0.9998747
[33,] 7.367098e-05 1.473420e-04 0.9999263
[34,] 6.251492e-05 1.250298e-04 0.9999375
[35,] 2.854213e-05 5.708427e-05 0.9999715
[36,] 1.903750e-05 3.807500e-05 0.9999810
[37,] 1.293286e-05 2.586572e-05 0.9999871
[38,] 2.253025e-05 4.506051e-05 0.9999775
[39,] 3.050633e-05 6.101265e-05 0.9999695
[40,] 3.857351e-03 7.714703e-03 0.9961426
[41,] 4.387755e-03 8.775510e-03 0.9956122
[42,] 2.998190e-02 5.996381e-02 0.9700181
> postscript(file="/var/www/html/rcomp/tmp/1eups1259332455.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/2vqhm1259332455.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/3uivh1259332455.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/4k8z21259332455.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/5pzni1259332455.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 = 57
Frequency = 1
1 2 3 4 5 6
-31.924957 -3.231841 -5.020852 6.282852 -49.331231 14.932635
7 8 9 10 11 12
-9.000839 13.299038 45.421090 25.871321 -37.118505 32.860196
13 14 15 16 17 18
-56.532469 -22.321475 18.055043 -29.967420 44.511277 -40.458695
19 20 21 22 23 24
17.584441 -43.869682 46.124234 -23.597488 52.074356 -9.808758
25 26 27 28 29 30
35.664922 13.279652 -30.631267 18.361978 1.126054 29.851573
31 32 33 34 35 36
29.631332 42.889840 24.415799 21.377072 17.351171 -37.822269
37 38 39 40 41 42
20.281846 38.054565 24.837142 -36.183059 -46.165297 3.016206
43 44 45 46 47 48
24.669903 13.299750 -54.166081 21.904733 -88.425535 17.724629
49 50 51 52 53 54
10.583212 78.780625 37.444892 -56.769873 78.488456 9.364381
55 56 57
-79.664918 -121.874117 -15.529587
> postscript(file="/var/www/html/rcomp/tmp/6wgmz1259332455.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 = 57
Frequency = 1
lag(myerror, k = 1) myerror
0 -31.924957 NA
1 -3.231841 -31.924957
2 -5.020852 -3.231841
3 6.282852 -5.020852
4 -49.331231 6.282852
5 14.932635 -49.331231
6 -9.000839 14.932635
7 13.299038 -9.000839
8 45.421090 13.299038
9 25.871321 45.421090
10 -37.118505 25.871321
11 32.860196 -37.118505
12 -56.532469 32.860196
13 -22.321475 -56.532469
14 18.055043 -22.321475
15 -29.967420 18.055043
16 44.511277 -29.967420
17 -40.458695 44.511277
18 17.584441 -40.458695
19 -43.869682 17.584441
20 46.124234 -43.869682
21 -23.597488 46.124234
22 52.074356 -23.597488
23 -9.808758 52.074356
24 35.664922 -9.808758
25 13.279652 35.664922
26 -30.631267 13.279652
27 18.361978 -30.631267
28 1.126054 18.361978
29 29.851573 1.126054
30 29.631332 29.851573
31 42.889840 29.631332
32 24.415799 42.889840
33 21.377072 24.415799
34 17.351171 21.377072
35 -37.822269 17.351171
36 20.281846 -37.822269
37 38.054565 20.281846
38 24.837142 38.054565
39 -36.183059 24.837142
40 -46.165297 -36.183059
41 3.016206 -46.165297
42 24.669903 3.016206
43 13.299750 24.669903
44 -54.166081 13.299750
45 21.904733 -54.166081
46 -88.425535 21.904733
47 17.724629 -88.425535
48 10.583212 17.724629
49 78.780625 10.583212
50 37.444892 78.780625
51 -56.769873 37.444892
52 78.488456 -56.769873
53 9.364381 78.488456
54 -79.664918 9.364381
55 -121.874117 -79.664918
56 -15.529587 -121.874117
57 NA -15.529587
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -3.231841 -31.924957
[2,] -5.020852 -3.231841
[3,] 6.282852 -5.020852
[4,] -49.331231 6.282852
[5,] 14.932635 -49.331231
[6,] -9.000839 14.932635
[7,] 13.299038 -9.000839
[8,] 45.421090 13.299038
[9,] 25.871321 45.421090
[10,] -37.118505 25.871321
[11,] 32.860196 -37.118505
[12,] -56.532469 32.860196
[13,] -22.321475 -56.532469
[14,] 18.055043 -22.321475
[15,] -29.967420 18.055043
[16,] 44.511277 -29.967420
[17,] -40.458695 44.511277
[18,] 17.584441 -40.458695
[19,] -43.869682 17.584441
[20,] 46.124234 -43.869682
[21,] -23.597488 46.124234
[22,] 52.074356 -23.597488
[23,] -9.808758 52.074356
[24,] 35.664922 -9.808758
[25,] 13.279652 35.664922
[26,] -30.631267 13.279652
[27,] 18.361978 -30.631267
[28,] 1.126054 18.361978
[29,] 29.851573 1.126054
[30,] 29.631332 29.851573
[31,] 42.889840 29.631332
[32,] 24.415799 42.889840
[33,] 21.377072 24.415799
[34,] 17.351171 21.377072
[35,] -37.822269 17.351171
[36,] 20.281846 -37.822269
[37,] 38.054565 20.281846
[38,] 24.837142 38.054565
[39,] -36.183059 24.837142
[40,] -46.165297 -36.183059
[41,] 3.016206 -46.165297
[42,] 24.669903 3.016206
[43,] 13.299750 24.669903
[44,] -54.166081 13.299750
[45,] 21.904733 -54.166081
[46,] -88.425535 21.904733
[47,] 17.724629 -88.425535
[48,] 10.583212 17.724629
[49,] 78.780625 10.583212
[50,] 37.444892 78.780625
[51,] -56.769873 37.444892
[52,] 78.488456 -56.769873
[53,] 9.364381 78.488456
[54,] -79.664918 9.364381
[55,] -121.874117 -79.664918
[56,] -15.529587 -121.874117
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -3.231841 -31.924957
2 -5.020852 -3.231841
3 6.282852 -5.020852
4 -49.331231 6.282852
5 14.932635 -49.331231
6 -9.000839 14.932635
7 13.299038 -9.000839
8 45.421090 13.299038
9 25.871321 45.421090
10 -37.118505 25.871321
11 32.860196 -37.118505
12 -56.532469 32.860196
13 -22.321475 -56.532469
14 18.055043 -22.321475
15 -29.967420 18.055043
16 44.511277 -29.967420
17 -40.458695 44.511277
18 17.584441 -40.458695
19 -43.869682 17.584441
20 46.124234 -43.869682
21 -23.597488 46.124234
22 52.074356 -23.597488
23 -9.808758 52.074356
24 35.664922 -9.808758
25 13.279652 35.664922
26 -30.631267 13.279652
27 18.361978 -30.631267
28 1.126054 18.361978
29 29.851573 1.126054
30 29.631332 29.851573
31 42.889840 29.631332
32 24.415799 42.889840
33 21.377072 24.415799
34 17.351171 21.377072
35 -37.822269 17.351171
36 20.281846 -37.822269
37 38.054565 20.281846
38 24.837142 38.054565
39 -36.183059 24.837142
40 -46.165297 -36.183059
41 3.016206 -46.165297
42 24.669903 3.016206
43 13.299750 24.669903
44 -54.166081 13.299750
45 21.904733 -54.166081
46 -88.425535 21.904733
47 17.724629 -88.425535
48 10.583212 17.724629
49 78.780625 10.583212
50 37.444892 78.780625
51 -56.769873 37.444892
52 78.488456 -56.769873
53 9.364381 78.488456
54 -79.664918 9.364381
55 -121.874117 -79.664918
56 -15.529587 -121.874117
> 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/7zjgc1259332455.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/8gfbp1259332455.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/9rk371259332455.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/10wfic1259332455.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/118m5a1259332455.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/12qo8c1259332455.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/13gqas1259332455.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/14e7321259332455.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/15xuku1259332455.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/16ohol1259332455.tab")
+ }
>
> system("convert tmp/1eups1259332455.ps tmp/1eups1259332455.png")
> system("convert tmp/2vqhm1259332455.ps tmp/2vqhm1259332455.png")
> system("convert tmp/3uivh1259332455.ps tmp/3uivh1259332455.png")
> system("convert tmp/4k8z21259332455.ps tmp/4k8z21259332455.png")
> system("convert tmp/5pzni1259332455.ps tmp/5pzni1259332455.png")
> system("convert tmp/6wgmz1259332455.ps tmp/6wgmz1259332455.png")
> system("convert tmp/7zjgc1259332455.ps tmp/7zjgc1259332455.png")
> system("convert tmp/8gfbp1259332455.ps tmp/8gfbp1259332455.png")
> system("convert tmp/9rk371259332455.ps tmp/9rk371259332455.png")
> system("convert tmp/10wfic1259332455.ps tmp/10wfic1259332455.png")
>
>
> proc.time()
user system elapsed
2.441 1.556 3.025