R version 2.15.2 (2012-10-26) -- "Trick or Treat"
Copyright (C) 2012 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i686-pc-linux-gnu (32-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> x <- array(list(655362,873127,1107897,1555964,1671159,1493308,2957796,2638691,1305669,1280496,921900,867888,652586,913831,1108544,1555827,1699283,1509458,3268975,2425016,1312703,1365498,934453,775019,651142,843192,1146766,1652601,1465906,1652734,2922334,2702805,1458956,1410363,1019279,936574,708917,885295,1099663,1576220,1487870,1488635,2882530,2677026,1404398,1344370,936865,872705,628151,953712,1160384,1400618,1661511,1495347,2918786,2775677,1407026,1370199,964526,850851,683118,847224,1073256,1514326,1503734,1507712,2865698,2788128,1391596,1366378,946295,859626),dim=c(1,72),dimnames=list(c('Y'),1:72))
> y <- array(NA,dim=c(1,72),dimnames=list(c('Y'),1:72))
> 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'
> library(lattice)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following object(s) are masked from 'package:base':
as.Date, 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 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 655362 1 0 0 0 0 0 0 0 0 0 0 1
2 873127 0 1 0 0 0 0 0 0 0 0 0 2
3 1107897 0 0 1 0 0 0 0 0 0 0 0 3
4 1555964 0 0 0 1 0 0 0 0 0 0 0 4
5 1671159 0 0 0 0 1 0 0 0 0 0 0 5
6 1493308 0 0 0 0 0 1 0 0 0 0 0 6
7 2957796 0 0 0 0 0 0 1 0 0 0 0 7
8 2638691 0 0 0 0 0 0 0 1 0 0 0 8
9 1305669 0 0 0 0 0 0 0 0 1 0 0 9
10 1280496 0 0 0 0 0 0 0 0 0 1 0 10
11 921900 0 0 0 0 0 0 0 0 0 0 1 11
12 867888 0 0 0 0 0 0 0 0 0 0 0 12
13 652586 1 0 0 0 0 0 0 0 0 0 0 13
14 913831 0 1 0 0 0 0 0 0 0 0 0 14
15 1108544 0 0 1 0 0 0 0 0 0 0 0 15
16 1555827 0 0 0 1 0 0 0 0 0 0 0 16
17 1699283 0 0 0 0 1 0 0 0 0 0 0 17
18 1509458 0 0 0 0 0 1 0 0 0 0 0 18
19 3268975 0 0 0 0 0 0 1 0 0 0 0 19
20 2425016 0 0 0 0 0 0 0 1 0 0 0 20
21 1312703 0 0 0 0 0 0 0 0 1 0 0 21
22 1365498 0 0 0 0 0 0 0 0 0 1 0 22
23 934453 0 0 0 0 0 0 0 0 0 0 1 23
24 775019 0 0 0 0 0 0 0 0 0 0 0 24
25 651142 1 0 0 0 0 0 0 0 0 0 0 25
26 843192 0 1 0 0 0 0 0 0 0 0 0 26
27 1146766 0 0 1 0 0 0 0 0 0 0 0 27
28 1652601 0 0 0 1 0 0 0 0 0 0 0 28
29 1465906 0 0 0 0 1 0 0 0 0 0 0 29
30 1652734 0 0 0 0 0 1 0 0 0 0 0 30
31 2922334 0 0 0 0 0 0 1 0 0 0 0 31
32 2702805 0 0 0 0 0 0 0 1 0 0 0 32
33 1458956 0 0 0 0 0 0 0 0 1 0 0 33
34 1410363 0 0 0 0 0 0 0 0 0 1 0 34
35 1019279 0 0 0 0 0 0 0 0 0 0 1 35
36 936574 0 0 0 0 0 0 0 0 0 0 0 36
37 708917 1 0 0 0 0 0 0 0 0 0 0 37
38 885295 0 1 0 0 0 0 0 0 0 0 0 38
39 1099663 0 0 1 0 0 0 0 0 0 0 0 39
40 1576220 0 0 0 1 0 0 0 0 0 0 0 40
41 1487870 0 0 0 0 1 0 0 0 0 0 0 41
42 1488635 0 0 0 0 0 1 0 0 0 0 0 42
43 2882530 0 0 0 0 0 0 1 0 0 0 0 43
44 2677026 0 0 0 0 0 0 0 1 0 0 0 44
45 1404398 0 0 0 0 0 0 0 0 1 0 0 45
46 1344370 0 0 0 0 0 0 0 0 0 1 0 46
47 936865 0 0 0 0 0 0 0 0 0 0 1 47
48 872705 0 0 0 0 0 0 0 0 0 0 0 48
49 628151 1 0 0 0 0 0 0 0 0 0 0 49
50 953712 0 1 0 0 0 0 0 0 0 0 0 50
51 1160384 0 0 1 0 0 0 0 0 0 0 0 51
52 1400618 0 0 0 1 0 0 0 0 0 0 0 52
53 1661511 0 0 0 0 1 0 0 0 0 0 0 53
54 1495347 0 0 0 0 0 1 0 0 0 0 0 54
55 2918786 0 0 0 0 0 0 1 0 0 0 0 55
56 2775677 0 0 0 0 0 0 0 1 0 0 0 56
57 1407026 0 0 0 0 0 0 0 0 1 0 0 57
58 1370199 0 0 0 0 0 0 0 0 0 1 0 58
59 964526 0 0 0 0 0 0 0 0 0 0 1 59
60 850851 0 0 0 0 0 0 0 0 0 0 0 60
61 683118 1 0 0 0 0 0 0 0 0 0 0 61
62 847224 0 1 0 0 0 0 0 0 0 0 0 62
63 1073256 0 0 1 0 0 0 0 0 0 0 0 63
64 1514326 0 0 0 1 0 0 0 0 0 0 0 64
65 1503734 0 0 0 0 1 0 0 0 0 0 0 65
66 1507712 0 0 0 0 0 1 0 0 0 0 0 66
67 2865698 0 0 0 0 0 0 1 0 0 0 0 67
68 2788128 0 0 0 0 0 0 0 1 0 0 0 68
69 1391596 0 0 0 0 0 0 0 0 1 0 0 69
70 1366378 0 0 0 0 0 0 0 0 0 1 0 70
71 946295 0 0 0 0 0 0 0 0 0 0 1 71
72 859626 0 0 0 0 0 0 0 0 0 0 0 72
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) M1 M2 M3 M4 M5
862212.2 -197694.3 25198.6 255262.2 681812.0 720838.6
M6 M7 M8 M9 M10 M11
663835.9 2108698.8 1807278.3 519487.9 495689.3 93400.4
t
-42.1
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-243632 -33513 -8567 28401 298864
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 862212.2 37852.3 22.778 < 2e-16 ***
M1 -197694.3 46348.0 -4.265 7.32e-05 ***
M2 25198.6 46300.3 0.544 0.5883
M3 255262.2 46257.1 5.518 8.02e-07 ***
M4 681812.0 46218.4 14.752 < 2e-16 ***
M5 720838.6 46184.2 15.608 < 2e-16 ***
M6 663835.9 46154.5 14.383 < 2e-16 ***
M7 2108698.8 46129.4 45.713 < 2e-16 ***
M8 1807278.2 46108.9 39.196 < 2e-16 ***
M9 519487.9 46092.9 11.270 2.41e-16 ***
M10 495689.3 46081.5 10.757 1.56e-15 ***
M11 93400.4 46074.6 2.027 0.0472 *
t -42.1 458.9 -0.092 0.9272
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 79800 on 59 degrees of freedom
Multiple R-squared: 0.9887, Adjusted R-squared: 0.9864
F-statistic: 430.1 on 12 and 59 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,] 0.0060830528 0.0121661057 0.9939169472
[2,] 0.0013110211 0.0026220422 0.9986889789
[3,] 0.0001775137 0.0003550274 0.9998224863
[4,] 0.6861641900 0.6276716199 0.3138358100
[5,] 0.9848276689 0.0303446622 0.0151723311
[6,] 0.9806457139 0.0387085721 0.0193542861
[7,] 0.9695836422 0.0608327155 0.0304163578
[8,] 0.9518245257 0.0963509485 0.0481754743
[9,] 0.9686169905 0.0627660190 0.0313830095
[10,] 0.9525732745 0.0948534510 0.0474267255
[11,] 0.9469393431 0.1061213137 0.0530606569
[12,] 0.9202702903 0.1594594195 0.0797297097
[13,] 0.9462596735 0.1074806530 0.0537403265
[14,] 0.9922792825 0.0154414350 0.0077207175
[15,] 0.9980272663 0.0039454674 0.0019727337
[16,] 0.9988947602 0.0022104795 0.0011052398
[17,] 0.9993332084 0.0013335833 0.0006667916
[18,] 0.9992632605 0.0014734789 0.0007367395
[19,] 0.9987591014 0.0024817972 0.0012408986
[20,] 0.9982933040 0.0034133921 0.0017066960
[21,] 0.9981090580 0.0037818841 0.0018909420
[22,] 0.9970278557 0.0059442887 0.0029721443
[23,] 0.9945212056 0.0109575888 0.0054787944
[24,] 0.9907750520 0.0184498960 0.0092249480
[25,] 0.9949156578 0.0101686844 0.0050843422
[26,] 0.9969153315 0.0061693369 0.0030846685
[27,] 0.9945827196 0.0108345609 0.0054172804
[28,] 0.9936489082 0.0127021837 0.0063510918
[29,] 0.9965318829 0.0069362342 0.0034681171
[30,] 0.9929557958 0.0140884083 0.0070442042
[31,] 0.9880149965 0.0239700070 0.0119850035
[32,] 0.9804010523 0.0391978953 0.0195989477
[33,] 0.9635693568 0.0728612865 0.0364306432
[34,] 0.9553755833 0.0892488334 0.0446244167
[35,] 0.9498357083 0.1003285834 0.0501642917
[36,] 0.9347898227 0.1304203546 0.0652101773
[37,] 0.9796129439 0.0407741122 0.0203870561
[38,] 0.9997634552 0.0004730896 0.0002365448
[39,] 0.9991284522 0.0017430956 0.0008715478
[40,] 0.9996527053 0.0006945895 0.0003472947
[41,] 0.9986267815 0.0027464369 0.0013732185
> postscript(file="/var/fisher/rcomp/tmp/13qte1355855454.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(x[,1], type='l', main='Actuals and Interpolation', ylab='value of Actuals and Interpolation (dots)', xlab='time or index')
> points(x[,1]-mysum$resid)
> grid()
> dev.off()
null device
1
> postscript(file="/var/fisher/rcomp/tmp/27bpe1355855454.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(mysum$resid, type='b', pch=19, main='Residuals', ylab='value of Residuals', xlab='time or index')
> grid()
> dev.off()
null device
1
> postscript(file="/var/fisher/rcomp/tmp/3wsam1355855454.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> hist(mysum$resid, main='Residual Histogram', xlab='values of Residuals')
> grid()
> dev.off()
null device
1
> postscript(file="/var/fisher/rcomp/tmp/42d7f1355855454.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> densityplot(~mysum$resid,col='black',main='Residual Density Plot', xlab='values of Residuals')
> dev.off()
null device
1
> postscript(file="/var/fisher/rcomp/tmp/50kvo1355855454.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> qqnorm(mysum$resid, main='Residual Normal Q-Q Plot')
> qqline(mysum$resid)
> grid()
> dev.off()
null device
1
> (myerror <- as.ts(mysum$resid))
Time Series:
Start = 1
End = 72
Frequency = 1
1 2 3 4 5 6
-9113.7798 -14199.6131 -9451.1131 12108.2202 88318.7202 -32487.4464
7 8 9 10 11 12
-12820.2798 -30462.6131 -75652.1131 -76984.4464 -33249.4464 6181.0536
13 14 15 16 17 18
-11384.5345 27009.6321 -8298.8679 12476.4655 116947.9655 -15832.2012
19 20 21 22 23 24
298863.9655 -243632.3679 -68112.8679 8522.7988 -20191.2012 -86182.7012
25 26 27 28 29 30
-12323.2893 -43124.1226 30428.3774 109755.7107 -115923.7893 127949.0440
31 32 33 34 35 36
-47271.7893 34661.8774 78645.3774 53893.0440 65140.0440 75877.5440
37 38 39 40 41 42
45956.9560 -515.8774 -16169.3774 33879.9560 -93454.5440 -35644.7107
43 44 45 46 47 48
-86570.5440 9388.1226 24592.6226 -11594.7107 -16768.7107 12513.7893
49 50 51 52 53 54
-34303.7988 68406.3679 45056.8679 -141216.7988 80691.7012 -28427.4655
55 56 57 58 59 60
-49809.2988 108544.3679 27725.8679 14739.5345 11397.5345 -8834.9655
61 62 63 64 65 66
21168.4464 -37576.3869 -41565.8869 -27003.5536 -76580.0536 -15557.2202
67 68 69 70 71 72
-102392.0536 121500.6131 12801.1131 11423.7798 -6328.2202 445.2798
> postscript(file="/var/fisher/rcomp/tmp/6ztq41355855454.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 72
Frequency = 1
lag(myerror, k = 1) myerror
0 -9113.7798 NA
1 -14199.6131 -9113.7798
2 -9451.1131 -14199.6131
3 12108.2202 -9451.1131
4 88318.7202 12108.2202
5 -32487.4464 88318.7202
6 -12820.2798 -32487.4464
7 -30462.6131 -12820.2798
8 -75652.1131 -30462.6131
9 -76984.4464 -75652.1131
10 -33249.4464 -76984.4464
11 6181.0536 -33249.4464
12 -11384.5345 6181.0536
13 27009.6321 -11384.5345
14 -8298.8679 27009.6321
15 12476.4655 -8298.8679
16 116947.9655 12476.4655
17 -15832.2012 116947.9655
18 298863.9655 -15832.2012
19 -243632.3679 298863.9655
20 -68112.8679 -243632.3679
21 8522.7988 -68112.8679
22 -20191.2012 8522.7988
23 -86182.7012 -20191.2012
24 -12323.2893 -86182.7012
25 -43124.1226 -12323.2893
26 30428.3774 -43124.1226
27 109755.7107 30428.3774
28 -115923.7893 109755.7107
29 127949.0440 -115923.7893
30 -47271.7893 127949.0440
31 34661.8774 -47271.7893
32 78645.3774 34661.8774
33 53893.0440 78645.3774
34 65140.0440 53893.0440
35 75877.5440 65140.0440
36 45956.9560 75877.5440
37 -515.8774 45956.9560
38 -16169.3774 -515.8774
39 33879.9560 -16169.3774
40 -93454.5440 33879.9560
41 -35644.7107 -93454.5440
42 -86570.5440 -35644.7107
43 9388.1226 -86570.5440
44 24592.6226 9388.1226
45 -11594.7107 24592.6226
46 -16768.7107 -11594.7107
47 12513.7893 -16768.7107
48 -34303.7988 12513.7893
49 68406.3679 -34303.7988
50 45056.8679 68406.3679
51 -141216.7988 45056.8679
52 80691.7012 -141216.7988
53 -28427.4655 80691.7012
54 -49809.2988 -28427.4655
55 108544.3679 -49809.2988
56 27725.8679 108544.3679
57 14739.5345 27725.8679
58 11397.5345 14739.5345
59 -8834.9655 11397.5345
60 21168.4464 -8834.9655
61 -37576.3869 21168.4464
62 -41565.8869 -37576.3869
63 -27003.5536 -41565.8869
64 -76580.0536 -27003.5536
65 -15557.2202 -76580.0536
66 -102392.0536 -15557.2202
67 121500.6131 -102392.0536
68 12801.1131 121500.6131
69 11423.7798 12801.1131
70 -6328.2202 11423.7798
71 445.2798 -6328.2202
72 NA 445.2798
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -14199.6131 -9113.7798
[2,] -9451.1131 -14199.6131
[3,] 12108.2202 -9451.1131
[4,] 88318.7202 12108.2202
[5,] -32487.4464 88318.7202
[6,] -12820.2798 -32487.4464
[7,] -30462.6131 -12820.2798
[8,] -75652.1131 -30462.6131
[9,] -76984.4464 -75652.1131
[10,] -33249.4464 -76984.4464
[11,] 6181.0536 -33249.4464
[12,] -11384.5345 6181.0536
[13,] 27009.6321 -11384.5345
[14,] -8298.8679 27009.6321
[15,] 12476.4655 -8298.8679
[16,] 116947.9655 12476.4655
[17,] -15832.2012 116947.9655
[18,] 298863.9655 -15832.2012
[19,] -243632.3679 298863.9655
[20,] -68112.8679 -243632.3679
[21,] 8522.7988 -68112.8679
[22,] -20191.2012 8522.7988
[23,] -86182.7012 -20191.2012
[24,] -12323.2893 -86182.7012
[25,] -43124.1226 -12323.2893
[26,] 30428.3774 -43124.1226
[27,] 109755.7107 30428.3774
[28,] -115923.7893 109755.7107
[29,] 127949.0440 -115923.7893
[30,] -47271.7893 127949.0440
[31,] 34661.8774 -47271.7893
[32,] 78645.3774 34661.8774
[33,] 53893.0440 78645.3774
[34,] 65140.0440 53893.0440
[35,] 75877.5440 65140.0440
[36,] 45956.9560 75877.5440
[37,] -515.8774 45956.9560
[38,] -16169.3774 -515.8774
[39,] 33879.9560 -16169.3774
[40,] -93454.5440 33879.9560
[41,] -35644.7107 -93454.5440
[42,] -86570.5440 -35644.7107
[43,] 9388.1226 -86570.5440
[44,] 24592.6226 9388.1226
[45,] -11594.7107 24592.6226
[46,] -16768.7107 -11594.7107
[47,] 12513.7893 -16768.7107
[48,] -34303.7988 12513.7893
[49,] 68406.3679 -34303.7988
[50,] 45056.8679 68406.3679
[51,] -141216.7988 45056.8679
[52,] 80691.7012 -141216.7988
[53,] -28427.4655 80691.7012
[54,] -49809.2988 -28427.4655
[55,] 108544.3679 -49809.2988
[56,] 27725.8679 108544.3679
[57,] 14739.5345 27725.8679
[58,] 11397.5345 14739.5345
[59,] -8834.9655 11397.5345
[60,] 21168.4464 -8834.9655
[61,] -37576.3869 21168.4464
[62,] -41565.8869 -37576.3869
[63,] -27003.5536 -41565.8869
[64,] -76580.0536 -27003.5536
[65,] -15557.2202 -76580.0536
[66,] -102392.0536 -15557.2202
[67,] 121500.6131 -102392.0536
[68,] 12801.1131 121500.6131
[69,] 11423.7798 12801.1131
[70,] -6328.2202 11423.7798
[71,] 445.2798 -6328.2202
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -14199.6131 -9113.7798
2 -9451.1131 -14199.6131
3 12108.2202 -9451.1131
4 88318.7202 12108.2202
5 -32487.4464 88318.7202
6 -12820.2798 -32487.4464
7 -30462.6131 -12820.2798
8 -75652.1131 -30462.6131
9 -76984.4464 -75652.1131
10 -33249.4464 -76984.4464
11 6181.0536 -33249.4464
12 -11384.5345 6181.0536
13 27009.6321 -11384.5345
14 -8298.8679 27009.6321
15 12476.4655 -8298.8679
16 116947.9655 12476.4655
17 -15832.2012 116947.9655
18 298863.9655 -15832.2012
19 -243632.3679 298863.9655
20 -68112.8679 -243632.3679
21 8522.7988 -68112.8679
22 -20191.2012 8522.7988
23 -86182.7012 -20191.2012
24 -12323.2893 -86182.7012
25 -43124.1226 -12323.2893
26 30428.3774 -43124.1226
27 109755.7107 30428.3774
28 -115923.7893 109755.7107
29 127949.0440 -115923.7893
30 -47271.7893 127949.0440
31 34661.8774 -47271.7893
32 78645.3774 34661.8774
33 53893.0440 78645.3774
34 65140.0440 53893.0440
35 75877.5440 65140.0440
36 45956.9560 75877.5440
37 -515.8774 45956.9560
38 -16169.3774 -515.8774
39 33879.9560 -16169.3774
40 -93454.5440 33879.9560
41 -35644.7107 -93454.5440
42 -86570.5440 -35644.7107
43 9388.1226 -86570.5440
44 24592.6226 9388.1226
45 -11594.7107 24592.6226
46 -16768.7107 -11594.7107
47 12513.7893 -16768.7107
48 -34303.7988 12513.7893
49 68406.3679 -34303.7988
50 45056.8679 68406.3679
51 -141216.7988 45056.8679
52 80691.7012 -141216.7988
53 -28427.4655 80691.7012
54 -49809.2988 -28427.4655
55 108544.3679 -49809.2988
56 27725.8679 108544.3679
57 14739.5345 27725.8679
58 11397.5345 14739.5345
59 -8834.9655 11397.5345
60 21168.4464 -8834.9655
61 -37576.3869 21168.4464
62 -41565.8869 -37576.3869
63 -27003.5536 -41565.8869
64 -76580.0536 -27003.5536
65 -15557.2202 -76580.0536
66 -102392.0536 -15557.2202
67 121500.6131 -102392.0536
68 12801.1131 121500.6131
69 11423.7798 12801.1131
70 -6328.2202 11423.7798
71 445.2798 -6328.2202
> 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/fisher/rcomp/tmp/7yhxu1355855454.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> acf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/fisher/rcomp/tmp/89r0c1355855454.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> pacf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Partial Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/fisher/rcomp/tmp/9v2wa1355855454.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
> plot(mylm, las = 1, sub='Residual Diagnostics')
> par(opar)
> dev.off()
null device
1
> if (n > n25) {
+ postscript(file="/var/fisher/rcomp/tmp/10274p1355855454.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
+ plot(kp3:nmkm3,gqarr[,2], main='Goldfeld-Quandt test',ylab='2-sided p-value',xlab='breakpoint')
+ grid()
+ dev.off()
+ }
null device
1
>
> #Note: the /var/fisher/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/fisher/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/fisher/rcomp/tmp/11l8qi1355855454.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/fisher/rcomp/tmp/12jpxz1355855454.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/fisher/rcomp/tmp/130h6y1355855454.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/fisher/rcomp/tmp/14nuiw1355855455.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/fisher/rcomp/tmp/1545hq1355855455.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/fisher/rcomp/tmp/16gc891355855455.tab")
+ }
>
> try(system("convert tmp/13qte1355855454.ps tmp/13qte1355855454.png",intern=TRUE))
character(0)
> try(system("convert tmp/27bpe1355855454.ps tmp/27bpe1355855454.png",intern=TRUE))
character(0)
> try(system("convert tmp/3wsam1355855454.ps tmp/3wsam1355855454.png",intern=TRUE))
character(0)
> try(system("convert tmp/42d7f1355855454.ps tmp/42d7f1355855454.png",intern=TRUE))
character(0)
> try(system("convert tmp/50kvo1355855454.ps tmp/50kvo1355855454.png",intern=TRUE))
character(0)
> try(system("convert tmp/6ztq41355855454.ps tmp/6ztq41355855454.png",intern=TRUE))
character(0)
> try(system("convert tmp/7yhxu1355855454.ps tmp/7yhxu1355855454.png",intern=TRUE))
character(0)
> try(system("convert tmp/89r0c1355855454.ps tmp/89r0c1355855454.png",intern=TRUE))
character(0)
> try(system("convert tmp/9v2wa1355855454.ps tmp/9v2wa1355855454.png",intern=TRUE))
character(0)
> try(system("convert tmp/10274p1355855454.ps tmp/10274p1355855454.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
6.474 1.724 8.246