R version 2.6.2 (2008-02-08)
Copyright (C) 2008 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
Natural language support but running in an English locale
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> x <- array(list(56421,53152,53536,52408,41454,38271,35306,26414,31917,38030,27534,18387,50556,43901,48572,43899,37532,40357,35489,29027,34485,42598,30306,26451,47460,50104,61465,53726,39477,43895,31481,29896,33842,39120,33702,25094,51442,45594,52518,48564,41745,49585,32747,33379,35645,37034,35681,20972,58552,54955,65540,51570,51145,46641,35704,33253,35193,41668,34865,21210,56126,49231,59723,48103,47472,50497,40059,34149,36860,46356,36577),dim=c(1,71),dimnames=list(c('x'),1:71))
> y <- array(NA,dim=c(1,71),dimnames=list(c('x'),1:71))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'Linear Trend'
> par2 = 'Include Monthly Dummies'
> par1 = '1'
> #'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
> #Author: Prof. Dr. P. Wessa
> #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #Technical description: Write here your technical program description (don't use hard returns!)
> library(lattice)
> library(lmtest)
Loading required package: zoo
> n25 <- 25 #minimum number of obs. for Goldfeld-Quandt test
> par1 <- as.numeric(par1)
> x <- t(y)
> k <- length(x[1,])
> n <- length(x[,1])
> x1 <- cbind(x[,par1], x[,1:k!=par1])
> mycolnames <- c(colnames(x)[par1], colnames(x)[1:k!=par1])
> colnames(x1) <- mycolnames #colnames(x)[par1]
> x <- x1
> if (par3 == 'First Differences'){
+ x2 <- array(0, dim=c(n-1,k), dimnames=list(1:(n-1), paste('(1-B)',colnames(x),sep='')))
+ for (i in 1:n-1) {
+ for (j in 1:k) {
+ x2[i,j] <- x[i+1,j] - x[i,j]
+ }
+ }
+ x <- x2
+ }
> if (par2 == 'Include Monthly Dummies'){
+ x2 <- array(0, dim=c(n,11), dimnames=list(1:n, paste('M', seq(1:11), sep ='')))
+ for (i in 1:11){
+ x2[seq(i,n,12),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> if (par2 == 'Include Quarterly Dummies'){
+ x2 <- array(0, dim=c(n,3), dimnames=list(1:n, paste('Q', seq(1:3), sep ='')))
+ for (i in 1:3){
+ x2[seq(i,n,4),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> k <- length(x[1,])
> if (par3 == 'Linear Trend'){
+ x <- cbind(x, c(1:n))
+ colnames(x)[k+1] <- 't'
+ }
> x
x M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 56421 1 0 0 0 0 0 0 0 0 0 0 1
2 53152 0 1 0 0 0 0 0 0 0 0 0 2
3 53536 0 0 1 0 0 0 0 0 0 0 0 3
4 52408 0 0 0 1 0 0 0 0 0 0 0 4
5 41454 0 0 0 0 1 0 0 0 0 0 0 5
6 38271 0 0 0 0 0 1 0 0 0 0 0 6
7 35306 0 0 0 0 0 0 1 0 0 0 0 7
8 26414 0 0 0 0 0 0 0 1 0 0 0 8
9 31917 0 0 0 0 0 0 0 0 1 0 0 9
10 38030 0 0 0 0 0 0 0 0 0 1 0 10
11 27534 0 0 0 0 0 0 0 0 0 0 1 11
12 18387 0 0 0 0 0 0 0 0 0 0 0 12
13 50556 1 0 0 0 0 0 0 0 0 0 0 13
14 43901 0 1 0 0 0 0 0 0 0 0 0 14
15 48572 0 0 1 0 0 0 0 0 0 0 0 15
16 43899 0 0 0 1 0 0 0 0 0 0 0 16
17 37532 0 0 0 0 1 0 0 0 0 0 0 17
18 40357 0 0 0 0 0 1 0 0 0 0 0 18
19 35489 0 0 0 0 0 0 1 0 0 0 0 19
20 29027 0 0 0 0 0 0 0 1 0 0 0 20
21 34485 0 0 0 0 0 0 0 0 1 0 0 21
22 42598 0 0 0 0 0 0 0 0 0 1 0 22
23 30306 0 0 0 0 0 0 0 0 0 0 1 23
24 26451 0 0 0 0 0 0 0 0 0 0 0 24
25 47460 1 0 0 0 0 0 0 0 0 0 0 25
26 50104 0 1 0 0 0 0 0 0 0 0 0 26
27 61465 0 0 1 0 0 0 0 0 0 0 0 27
28 53726 0 0 0 1 0 0 0 0 0 0 0 28
29 39477 0 0 0 0 1 0 0 0 0 0 0 29
30 43895 0 0 0 0 0 1 0 0 0 0 0 30
31 31481 0 0 0 0 0 0 1 0 0 0 0 31
32 29896 0 0 0 0 0 0 0 1 0 0 0 32
33 33842 0 0 0 0 0 0 0 0 1 0 0 33
34 39120 0 0 0 0 0 0 0 0 0 1 0 34
35 33702 0 0 0 0 0 0 0 0 0 0 1 35
36 25094 0 0 0 0 0 0 0 0 0 0 0 36
37 51442 1 0 0 0 0 0 0 0 0 0 0 37
38 45594 0 1 0 0 0 0 0 0 0 0 0 38
39 52518 0 0 1 0 0 0 0 0 0 0 0 39
40 48564 0 0 0 1 0 0 0 0 0 0 0 40
41 41745 0 0 0 0 1 0 0 0 0 0 0 41
42 49585 0 0 0 0 0 1 0 0 0 0 0 42
43 32747 0 0 0 0 0 0 1 0 0 0 0 43
44 33379 0 0 0 0 0 0 0 1 0 0 0 44
45 35645 0 0 0 0 0 0 0 0 1 0 0 45
46 37034 0 0 0 0 0 0 0 0 0 1 0 46
47 35681 0 0 0 0 0 0 0 0 0 0 1 47
48 20972 0 0 0 0 0 0 0 0 0 0 0 48
49 58552 1 0 0 0 0 0 0 0 0 0 0 49
50 54955 0 1 0 0 0 0 0 0 0 0 0 50
51 65540 0 0 1 0 0 0 0 0 0 0 0 51
52 51570 0 0 0 1 0 0 0 0 0 0 0 52
53 51145 0 0 0 0 1 0 0 0 0 0 0 53
54 46641 0 0 0 0 0 1 0 0 0 0 0 54
55 35704 0 0 0 0 0 0 1 0 0 0 0 55
56 33253 0 0 0 0 0 0 0 1 0 0 0 56
57 35193 0 0 0 0 0 0 0 0 1 0 0 57
58 41668 0 0 0 0 0 0 0 0 0 1 0 58
59 34865 0 0 0 0 0 0 0 0 0 0 1 59
60 21210 0 0 0 0 0 0 0 0 0 0 0 60
61 56126 1 0 0 0 0 0 0 0 0 0 0 61
62 49231 0 1 0 0 0 0 0 0 0 0 0 62
63 59723 0 0 1 0 0 0 0 0 0 0 0 63
64 48103 0 0 0 1 0 0 0 0 0 0 0 64
65 47472 0 0 0 0 1 0 0 0 0 0 0 65
66 50497 0 0 0 0 0 1 0 0 0 0 0 66
67 40059 0 0 0 0 0 0 1 0 0 0 0 67
68 34149 0 0 0 0 0 0 0 1 0 0 0 68
69 36860 0 0 0 0 0 0 0 0 1 0 0 69
70 46356 0 0 0 0 0 0 0 0 0 1 0 70
71 36577 0 0 0 0 0 0 0 0 0 0 1 71
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) M1 M2 M3 M4 M5
18935.82 31487.67 27454.14 34760.11 27482.59 20811.56
M6 M7 M8 M9 M10 M11
22451.53 12611.34 8403.15 11943.62 17990.76 10203.73
t
96.86
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-6577 -2589 -206 2006 6904
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 18935.82 1714.89 11.042 6.94e-16 ***
M1 31487.67 2102.76 14.974 < 2e-16 ***
M2 27454.14 2101.87 13.062 < 2e-16 ***
M3 34760.11 2101.19 16.543 < 2e-16 ***
M4 27482.59 2100.69 13.083 < 2e-16 ***
M5 20811.56 2100.40 9.908 4.36e-14 ***
M6 22451.53 2100.30 10.690 2.47e-15 ***
M7 12611.34 2100.40 6.004 1.34e-07 ***
M8 8403.15 2100.69 4.000 0.000182 ***
M9 11943.62 2101.19 5.684 4.49e-07 ***
M10 17990.76 2101.87 8.559 7.12e-12 ***
M11 10203.73 2102.76 4.853 9.56e-06 ***
t 96.86 20.31 4.769 1.29e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 3469 on 58 degrees of freedom
Multiple R-squared: 0.9075, Adjusted R-squared: 0.8884
F-statistic: 47.44 on 12 and 58 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.10367782 0.20735564 0.896322179
[2,] 0.09219828 0.18439657 0.907801717
[3,] 0.44205844 0.88411688 0.557941558
[4,] 0.46092518 0.92185036 0.539074821
[5,] 0.54935140 0.90129721 0.450648605
[6,] 0.57409917 0.85180167 0.425900833
[7,] 0.66931623 0.66136754 0.330683771
[8,] 0.64102059 0.71795883 0.358979413
[9,] 0.81829002 0.36341995 0.181709977
[10,] 0.82996087 0.34007826 0.170039131
[11,] 0.80287738 0.39424524 0.197122622
[12,] 0.94278962 0.11442076 0.057210382
[13,] 0.96778297 0.06443406 0.032217029
[14,] 0.95908358 0.08183283 0.040916417
[15,] 0.94708540 0.10582920 0.052914599
[16,] 0.93810288 0.12379423 0.061897115
[17,] 0.90953276 0.18093448 0.090467241
[18,] 0.87027131 0.25945737 0.129728686
[19,] 0.82425633 0.35148733 0.175743667
[20,] 0.78879083 0.42241834 0.211209170
[21,] 0.81382912 0.37234175 0.186170875
[22,] 0.77718451 0.44563099 0.222815494
[23,] 0.78083297 0.43833406 0.219167030
[24,] 0.88047385 0.23905230 0.119526152
[25,] 0.83449432 0.33101135 0.165505676
[26,] 0.88031178 0.23937645 0.119688225
[27,] 0.89682427 0.20635147 0.103175733
[28,] 0.89495547 0.21008905 0.105044525
[29,] 0.85803773 0.28392453 0.141962265
[30,] 0.79776154 0.40447691 0.202238457
[31,] 0.89322031 0.21355938 0.106779689
[32,] 0.84991202 0.30017595 0.150087976
[33,] 0.78985162 0.42029676 0.210148379
[34,] 0.74910088 0.50179825 0.250899125
[35,] 0.80394035 0.39211931 0.196059653
[36,] 0.91746753 0.16506494 0.082532468
[37,] 0.93912640 0.12174719 0.060873595
[38,] 0.99410503 0.01178993 0.005894966
[39,] 0.98190592 0.03618816 0.018094078
[40,] 0.96212200 0.07575599 0.037877997
> postscript(file="/var/www/html/rcomp/tmp/1pzuq1210175771.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/2w9521210175771.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/3p8su1210175771.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/4ya431210175771.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/5d2hv1210175771.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 = 71
Frequency = 1
1 2 3 4 5 6
5900.64815 6568.31481 -450.51852 5602.14815 1222.31481 -3697.51852
7 8 9 10 11 12
3080.81481 -1699.85185 165.81481 134.81481 -2671.01852 -1711.14815
13 14 15 16 17 18
-1126.67778 -3845.01111 -6576.84444 -4069.17778 -3862.01111 -2773.84444
19 20 21 22 23 24
2101.48889 -249.17778 1571.48889 3540.48889 -1061.34444 5190.52593
25 26 27 28 29 30
-5385.00370 1195.66296 5153.82963 4595.49630 -3079.33704 -398.17037
31 32 33 34 35 36
-3068.83704 -542.50370 -233.83704 -1099.83704 1172.32963 2671.20000
37 38 39 40 41 42
-2565.32963 -4476.66296 -4955.49630 -1728.82963 -1973.66296 4129.50370
43 44 45 46 47 48
-2965.16296 1778.17037 406.83704 -4348.16296 1989.00370 -2613.12593
49 50 51 52 53 54
3382.34444 3722.01111 6904.17778 114.84444 6264.01111 23.17778
55 56 57 58 59 60
-1170.48889 489.84444 -1207.48889 -876.48889 10.67778 -3537.45185
61 62 63 64 65 66
-205.98148 -3164.31481 -75.14815 -4514.48148 1428.68519 2716.85185
67 68 69 70 71
2022.18519 223.51852 -702.81481 2649.18519 560.35185
> postscript(file="/var/www/html/rcomp/tmp/6c0111210175771.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 = 71
Frequency = 1
lag(myerror, k = 1) myerror
0 5900.64815 NA
1 6568.31481 5900.64815
2 -450.51852 6568.31481
3 5602.14815 -450.51852
4 1222.31481 5602.14815
5 -3697.51852 1222.31481
6 3080.81481 -3697.51852
7 -1699.85185 3080.81481
8 165.81481 -1699.85185
9 134.81481 165.81481
10 -2671.01852 134.81481
11 -1711.14815 -2671.01852
12 -1126.67778 -1711.14815
13 -3845.01111 -1126.67778
14 -6576.84444 -3845.01111
15 -4069.17778 -6576.84444
16 -3862.01111 -4069.17778
17 -2773.84444 -3862.01111
18 2101.48889 -2773.84444
19 -249.17778 2101.48889
20 1571.48889 -249.17778
21 3540.48889 1571.48889
22 -1061.34444 3540.48889
23 5190.52593 -1061.34444
24 -5385.00370 5190.52593
25 1195.66296 -5385.00370
26 5153.82963 1195.66296
27 4595.49630 5153.82963
28 -3079.33704 4595.49630
29 -398.17037 -3079.33704
30 -3068.83704 -398.17037
31 -542.50370 -3068.83704
32 -233.83704 -542.50370
33 -1099.83704 -233.83704
34 1172.32963 -1099.83704
35 2671.20000 1172.32963
36 -2565.32963 2671.20000
37 -4476.66296 -2565.32963
38 -4955.49630 -4476.66296
39 -1728.82963 -4955.49630
40 -1973.66296 -1728.82963
41 4129.50370 -1973.66296
42 -2965.16296 4129.50370
43 1778.17037 -2965.16296
44 406.83704 1778.17037
45 -4348.16296 406.83704
46 1989.00370 -4348.16296
47 -2613.12593 1989.00370
48 3382.34444 -2613.12593
49 3722.01111 3382.34444
50 6904.17778 3722.01111
51 114.84444 6904.17778
52 6264.01111 114.84444
53 23.17778 6264.01111
54 -1170.48889 23.17778
55 489.84444 -1170.48889
56 -1207.48889 489.84444
57 -876.48889 -1207.48889
58 10.67778 -876.48889
59 -3537.45185 10.67778
60 -205.98148 -3537.45185
61 -3164.31481 -205.98148
62 -75.14815 -3164.31481
63 -4514.48148 -75.14815
64 1428.68519 -4514.48148
65 2716.85185 1428.68519
66 2022.18519 2716.85185
67 223.51852 2022.18519
68 -702.81481 223.51852
69 2649.18519 -702.81481
70 560.35185 2649.18519
71 NA 560.35185
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 6568.31481 5900.64815
[2,] -450.51852 6568.31481
[3,] 5602.14815 -450.51852
[4,] 1222.31481 5602.14815
[5,] -3697.51852 1222.31481
[6,] 3080.81481 -3697.51852
[7,] -1699.85185 3080.81481
[8,] 165.81481 -1699.85185
[9,] 134.81481 165.81481
[10,] -2671.01852 134.81481
[11,] -1711.14815 -2671.01852
[12,] -1126.67778 -1711.14815
[13,] -3845.01111 -1126.67778
[14,] -6576.84444 -3845.01111
[15,] -4069.17778 -6576.84444
[16,] -3862.01111 -4069.17778
[17,] -2773.84444 -3862.01111
[18,] 2101.48889 -2773.84444
[19,] -249.17778 2101.48889
[20,] 1571.48889 -249.17778
[21,] 3540.48889 1571.48889
[22,] -1061.34444 3540.48889
[23,] 5190.52593 -1061.34444
[24,] -5385.00370 5190.52593
[25,] 1195.66296 -5385.00370
[26,] 5153.82963 1195.66296
[27,] 4595.49630 5153.82963
[28,] -3079.33704 4595.49630
[29,] -398.17037 -3079.33704
[30,] -3068.83704 -398.17037
[31,] -542.50370 -3068.83704
[32,] -233.83704 -542.50370
[33,] -1099.83704 -233.83704
[34,] 1172.32963 -1099.83704
[35,] 2671.20000 1172.32963
[36,] -2565.32963 2671.20000
[37,] -4476.66296 -2565.32963
[38,] -4955.49630 -4476.66296
[39,] -1728.82963 -4955.49630
[40,] -1973.66296 -1728.82963
[41,] 4129.50370 -1973.66296
[42,] -2965.16296 4129.50370
[43,] 1778.17037 -2965.16296
[44,] 406.83704 1778.17037
[45,] -4348.16296 406.83704
[46,] 1989.00370 -4348.16296
[47,] -2613.12593 1989.00370
[48,] 3382.34444 -2613.12593
[49,] 3722.01111 3382.34444
[50,] 6904.17778 3722.01111
[51,] 114.84444 6904.17778
[52,] 6264.01111 114.84444
[53,] 23.17778 6264.01111
[54,] -1170.48889 23.17778
[55,] 489.84444 -1170.48889
[56,] -1207.48889 489.84444
[57,] -876.48889 -1207.48889
[58,] 10.67778 -876.48889
[59,] -3537.45185 10.67778
[60,] -205.98148 -3537.45185
[61,] -3164.31481 -205.98148
[62,] -75.14815 -3164.31481
[63,] -4514.48148 -75.14815
[64,] 1428.68519 -4514.48148
[65,] 2716.85185 1428.68519
[66,] 2022.18519 2716.85185
[67,] 223.51852 2022.18519
[68,] -702.81481 223.51852
[69,] 2649.18519 -702.81481
[70,] 560.35185 2649.18519
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 6568.31481 5900.64815
2 -450.51852 6568.31481
3 5602.14815 -450.51852
4 1222.31481 5602.14815
5 -3697.51852 1222.31481
6 3080.81481 -3697.51852
7 -1699.85185 3080.81481
8 165.81481 -1699.85185
9 134.81481 165.81481
10 -2671.01852 134.81481
11 -1711.14815 -2671.01852
12 -1126.67778 -1711.14815
13 -3845.01111 -1126.67778
14 -6576.84444 -3845.01111
15 -4069.17778 -6576.84444
16 -3862.01111 -4069.17778
17 -2773.84444 -3862.01111
18 2101.48889 -2773.84444
19 -249.17778 2101.48889
20 1571.48889 -249.17778
21 3540.48889 1571.48889
22 -1061.34444 3540.48889
23 5190.52593 -1061.34444
24 -5385.00370 5190.52593
25 1195.66296 -5385.00370
26 5153.82963 1195.66296
27 4595.49630 5153.82963
28 -3079.33704 4595.49630
29 -398.17037 -3079.33704
30 -3068.83704 -398.17037
31 -542.50370 -3068.83704
32 -233.83704 -542.50370
33 -1099.83704 -233.83704
34 1172.32963 -1099.83704
35 2671.20000 1172.32963
36 -2565.32963 2671.20000
37 -4476.66296 -2565.32963
38 -4955.49630 -4476.66296
39 -1728.82963 -4955.49630
40 -1973.66296 -1728.82963
41 4129.50370 -1973.66296
42 -2965.16296 4129.50370
43 1778.17037 -2965.16296
44 406.83704 1778.17037
45 -4348.16296 406.83704
46 1989.00370 -4348.16296
47 -2613.12593 1989.00370
48 3382.34444 -2613.12593
49 3722.01111 3382.34444
50 6904.17778 3722.01111
51 114.84444 6904.17778
52 6264.01111 114.84444
53 23.17778 6264.01111
54 -1170.48889 23.17778
55 489.84444 -1170.48889
56 -1207.48889 489.84444
57 -876.48889 -1207.48889
58 10.67778 -876.48889
59 -3537.45185 10.67778
60 -205.98148 -3537.45185
61 -3164.31481 -205.98148
62 -75.14815 -3164.31481
63 -4514.48148 -75.14815
64 1428.68519 -4514.48148
65 2716.85185 1428.68519
66 2022.18519 2716.85185
67 223.51852 2022.18519
68 -702.81481 223.51852
69 2649.18519 -702.81481
70 560.35185 2649.18519
> 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/7rese1210175771.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/8cq411210175771.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/980vf1210175771.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/10vgh21210175771.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/11wajj1210175771.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/12ozyi1210175771.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/13v87f1210175771.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/14kdz51210175771.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/15m4hr1210175772.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/1698vb1210175772.tab")
+ }
>
> system("convert tmp/1pzuq1210175771.ps tmp/1pzuq1210175771.png")
> system("convert tmp/2w9521210175771.ps tmp/2w9521210175771.png")
> system("convert tmp/3p8su1210175771.ps tmp/3p8su1210175771.png")
> system("convert tmp/4ya431210175771.ps tmp/4ya431210175771.png")
> system("convert tmp/5d2hv1210175771.ps tmp/5d2hv1210175771.png")
> system("convert tmp/6c0111210175771.ps tmp/6c0111210175771.png")
> system("convert tmp/7rese1210175771.ps tmp/7rese1210175771.png")
> system("convert tmp/8cq411210175771.ps tmp/8cq411210175771.png")
> system("convert tmp/980vf1210175771.ps tmp/980vf1210175771.png")
> system("convert tmp/10vgh21210175771.ps tmp/10vgh21210175771.png")
>
>
> proc.time()
user system elapsed
5.466 2.804 5.831