# R SCRIPT FOR STATISTICAL COMPUTING ANALYSIS # # # U.S. MONEY SUPPLY, PRICE LEVEL AND GDP DATA # library(lmtest) #usqdata <- read.table("e:\\DSLMHTML\\STATCOMP\\statcomp.tab",header=TRUE) usqdata <- read.table("statcomp.tab",header=TRUE) names(usqdata) attach(usqdata) # # Set up variables in the object usgdata as time-series USGDP <- ts(usqdata$USGDP,start=c(1959,1),end=c(2010,1),frequency=4) US3MTBR <- ts(usqdata$US3MTBR,start=c(1959,1),end=c(2010,1),frequency=4) USIPD <- ts(usqdata$USIPD,start=c(1959,1),end=c(2010,1),frequency=4) USCPI <- ts(usqdata$USCPI,start=c(1959,1),end=c(2010,1),frequency=4) USM1 <- ts(usqdata$USM1,start=c(1959,1),end=c(2010,1),frequency=4) USM2 <- ts(usqdata$USM2,start=c(1959,1),end=c(2010,1),frequency=4) # USCPIL4 <- lag(USCPI, -4) USIPDL4 <- lag(USIPD, -4) YYGCPI <- 100*(USCPI - USCPIL4)/USCPIL4 YYGIPD <- 100*(USIPD - USIPDL4)/USIPDL4 # plot(USCPI,type="l", xlab="Quarterly Data", ylab="Index -- 2005 = 100", main="United States Price Level") lines(USIPD,lty=2) legend(1960,110,c("CPI","GDP Deflator"),lty=c(1,2)) # plot(YYGCPI,type="l", xlab="Quarterly Data", ylab="Year-Over-Year Percentage Change", main="United States Year-Over-Year Inflation") lines(YYGIPD,lty=2) legend(1990,14,c("CPI","GDP Deflator"),lty=c(1,2)) # # Construct Kernel Densities for the two inflation rate # series and then plot them # DENYYGCPI <- density(YYGCPI) DENYYGIPD <- density(YYGIPD) # plot(DENYYGCPI, ylim=c(0,.25), main="Kernel Density Estimates of U.S. Inflation", sub="Year-over-Year") lines(DENYYGIPD,lty=2) legend(7,.22,c("CPI","GDP Deflator"),lty=c(1,2)) # # Examine details of the two inflation rate series # summary(YYGCPI) # Functions to give details of YYGCPI mean(YYGCPI) median(YYGCPI) var(YYGCPI) sd(YYGCPI) mean(YYGIPD) median(YYGIPD) var(YYGIPD) sd(YYGIPD) quantile(YYGCPI, .25) quantile(YYGCPI, .75) quantile(YYGCPI, .75) - quantile(YYGCPI, .25) # inter-quartile range max(YYGCPI) min(YYGCPI) max(YYGCPI) - min(YYGCPI) range(YYGCPI) # summary(YYGIPD) # Details of YYGIPD quantile(YYGIPD, .75) - quantile(YYGIPD, .25) # inter-quartile range range(YYGIPD) # # Estimate Demand Function for M1 in Logarithms # USRGDP <- 100*USGDP/USIPD # Calculate real GDP USRM1 <- 100*USM1/USIPD # Calculate real M1 LUSRGDP <- log(USRGDP) # Natural logarithm of real GDP LUSRM1 <- log(USRM1) # Natural logarithm of real M1 # dmreg <- lm(LUSRM1 ~ US3MTBR + LUSRGDP) # Demand for Money Regression -- 1959:Q1 through 2010:Q1 summary(dmreg) # Present basic regression results dmoutput <- summary(dmreg) SSR <- deviance(dmreg) DMCOEF <- dmreg$coefficients # regression coeficients DF <- dmreg$df DMFIT <- dmreg$fitted.values DMRES <- dmreg$residuals s <- dmoutput$sigma # standard error of regression residuals dmRSQ <- dmoutput$r.squared dmCovMat <- s^2*dmoutput$cov # variance-covariance matrix of coefficients } # DMRES <- ts(DMRES,start=c(1959,1),end=c(2010,1),frequency=4) # set up as DMFIT <- ts(DMFIT,start=c(1959,1),end=c(2010,1),frequency=4) # time series # plot(LUSRM1,type="l", xlab="Quarterly Data", ylab="Logarithms", main="Actual and Fitted: Demand for Money Regression") lines(DMFIT,lty=2) legend(1960,7.3,c("Actual","fitted"),lty=c(1,2)) # plot(DMRES,type="l", xlab="Quarterly Data", ylab=" ", main="Demand for Money Regression Residuals") # # Perform the Breusch-Pagan test for heteroskedasticity of the residuals # bptest(dmreg) bptest(dmreg,studentize=TRUE) bptest(dmreg,studentize=FALSE) # # Perform tests provided by the R-statistical program for serial correlation # in the residuals of the original demand-for-money regression, specifying # 4 lags of those residuals to be used in performing the tests # bgtest(dmreg,order=4) bgtest(dmreg,order=4,type="Chisq") bgtest(dmreg,order=4,type="F") Box.test(DMRES, lag = 4, type = "Ljung") # # Do the LM based test specified by Maddala using the Davidson and MacKinnon # view that unavailable lagged values should be set to zero # # Take lag of the residuals of the above regression # DMRESL1 <- lag(DMRES, -1) DMRESL2 <- lag(DMRES, -2) DMRESL3 <- lag(DMRES, -3) DMRESL4 <- lag(DMRES, -4) # # Put all time-series in a data-frame (i.e. worksheet form) which is ABSOLUTELY NECESSARY # to ensure that the lagged values are properly organized in relation to each other and # the relevant unlagged series. # BASEDAT <- ts.union(LUSRM1,US3MTBR,LUSRGDP,DMRES,DMRESL1,DMRESL2,DMRESL3,DMRESL4,USCPI, USIPD,YYGCPI,YYGIPD) # This worksheet can be examined using the command # edit(BASEDAT) # # Follow advice of Davidson and MacKinnon and replace missing lagged values of the # residuals with the number 0 --- the elements of BASEDAT are accessed by adding # the code fragment [row,column] to the name, where the columns are given by # the ordering of the variables in the ts.union function above. # BASEDAT[1,5] <- 0 #Assign the NA observations a value of 0 BASEDAT[1,6] <- 0 BASEDAT[1,7] <- 0 BASEDAT[1,8] <- 0 BASEDAT[2,6] <- 0 BASEDAT[2,7] <- 0 BASEDAT[2,8] <- 0 BASEDAT[3,7] <- 0 BASEDAT[3,8] <- 0 BASEDAT[4,8] <- 0 # # Now regress the residuals of the demand for money regression on the independent variables # in that regression plus the above constructed lags of those residuals. It is essential # to pass to the function the argument data=BASEMAT so it will use the variables in a # properly organized way. # uresreg <- lm(DMRES ~ DMRESL1 + DMRESL2 + DMRESL3 + DMRESL4 + US3MTBR + LUSRGDP, data=BASEDAT) summary(uresreg) # unrestricted residuals-regression uresregres <- uresreg$residuals # extract the regression residuals uresregresq <- uresregres^2 # square those residuals SSEU <- sum(uresregresq) # calculate the sum of squared errors of this unrestricted # regression # # Now run the above regression with the lagged residuals excluded -- this is the restricted rresreg <- lm(DMRES ~ US3MTBR + LUSRGDP, data=BASEDAT) # regression summary(rresreg) # restricted residuals-regression rresregres <- rresreg$residuals # residuals of above restricted regression rresregresq <- rresregres^2 # square those residuals SSER <- sum(rresregresq) # sum of squared errors of this restricted regression # URRSQ <- summary(uresreg)$r.squared # extract RSQ of unrestricted regression URRSQ # and print them RRRSQ <- summary(rresreg)$r.squared # extract RSQ of restricted regression RRRSQ # and print them # FNUM <- (SSER - SSEU)/4 # calculate the numerator of F-statistic SSER SSEU DMRESSQ <- DMRES^2 # calculate the sum of squares of residuals RRTSS <- sum(DMRESSQ) # of the original demand-for-money Regression RRTSS # FNUM # print the numerator of the F-Statistic FDEN <- SSEU/uresreg$df # calculate the denominator of F-statistic FDEN # print the denominator of the F-statistic FSTAT <- FNUM/FDEN # calculate the F-Statistic FSTAT # print the F-Statistic # qf(.99999, 4, 198) # F-Value with .00001 right-tail probability # LMCHISQSTAT <- 4*FSTAT # calculate the Chi-square statistic based on the F-Statistic LMCHISQSTAT # print out this Chi-square statistic # BGCHISQ <- 205*URRSQ # Calculate the Breusch Godfrey Chi-square statistic based on BGCHISQ # the R-Square of the unrestricted regression # qchisq(.99999, 4) # Chi-Square Value with .00001 right-tail probabiity # pf(10, 4, 198) pchisq(10, 4) # # BOOTSTRAP REGRESSION COEFFICIENTS # dmreg <- lm(LUSRM1 ~ US3MTBR + LUSRGDP) summary(dmreg) ORIGFITTED <- dmreg$fitted.values ORIGRESID <- dmreg$residuals # COEFSCON <- c(1:1000)-c(1:1000) COEFSTBR <- c(1:1000)-c(1:1000) COEFSGDP <- c(1:1000)-c(1:1000) # for (i in 1:1000) { NEWRESID <- sample(ORIGRESID) NEWLRM1 <- ORIGFITTED + NEWRESID newdmreg <- lm(NEWLRM1 ~ US3MTBR + LUSRGDP) NEWCOEFS <- newdmreg$coefficients COEFSCON[i] <- NEWCOEFS[1] COEFSTBR[i] <- NEWCOEFS[2] COEFSGDP[i] <- NEWCOEFS[3] } # summary(COEFSCON) summary(COEFSTBR) summary(COEFSGDP) # DENCOEFSCON <- density(COEFSCON) DENCOEFSTBR <- density(COEFSTBR) DENCOEFSGDP <- density(COEFSGDP) # plot(DENCOEFSCON, main="Kernel Density: Constant Term Coefficient") # plot(DENCOEFSTBR, main="Kernel Density: TBill Rate Coefficient") # plot(DENCOEFSGDP, main="Kernel Density: Log Real GDP Coefficient") #