Zelig September 2009

zelig@lists.gking.harvard.edu

10 participants
11 discussions

by Donald Braman

Dear Zelig Folk, For a few applications, I'd like to run a single simulation. I thought I might do this by setting num=1, but if I set num to anything less than 8, I get the following error: > Error in qi.polr(object[[1]], simpar = simpar, x = as.matrix(x), x1 = if > (!is.null(x1)) as.matrix(x1)) : subscript out of bounds Except when I set num to either 1 or 2, in which case I get the following error: Error in matrix(rnorm(p * n), n) : invalid 'ncol' value (too large or NA) It's not all that important to what I'm doing, but I'm curious if anyone knows why this happens or if there's anything that can be done about it. Cheers, Don

14 years, 6 months

quick question about error message

by Prashant

Hi, I'm trying to run matchit followed by zelig to calculate the ATT: After running z.out <- zelig(@@@), I get the following error message: Error in `[.data.frame`(d, , all.vars(as.expression(formula))) : undefined columns selected I'm not sure why I'm getting this error message--afterwards, I want to run: x.out1 <- setx(z.out1, data = match.data(m.outf, "treat"), cond = TRUE) s.out1 <- sim(z.out1, x = x.out1) Any help at all would be most appreciated! best, Prashant P.S. Specifically, I had run the following code: ps <- glm(Tr ~ female + urban + father_edu + mother_edu + siblings + hs_type_best + kuaiban + wangjie + sx_non_mil12 + I(age^2) + I(urban*age) + I(urban*mother_edu) + I(urban*sx_non_mil12) + I(female*sx_non_mil12) + I(urban*wangjie) + I(female*wangjie), family=binomial, data=my.data) X = cbind(age, female, urban, father_edu, mother_edu, siblings, hs_type_best, kuaiban, wangjie, sx_non_mil12, I(age^2), I(urban*age), I(urban*mother_edu), I(urban*sx_non_mil12), I(female*sx_non_mil12), I(urban*wangjie), I(female*wangjie), ps$fitted) m.outf <- matchit(Tr ~ X, data = rdmatchAC.lk12.range10, method = "genetic", ratio = 2, replace = TRUE, pop.size=1000, discard = "control") z.out1 <- zelig(Y1 ~ Tr + X + gk_score, model = "ls", data = match.data(m.outf, "control")) x.out1 <- setx(z.out1, data = match.data(m.outf, "treat"), cond = TRUE) s.out1 <- sim(z.out1, x = x.out1)

14 years, 6 months

Mixed effects regression model

by Sara Dias

Hey all, In my research I need to adjust a gamma mixed regression model. I'm using Zelig package. z<-zelig(y~a+d1+d2+d3+tag(1|CODE), data=los0, model="gamma.mixed") But I'm having the following problem: In mer_finalize(ans) : singular convergence (7) Can anyone explain me what's happening? What can I do, to do not have this error? Thank you very much. Sara Dias

14 years, 7 months

help with variable recoding in R

by Casey Klofstad

Hey all- I'm sure this is easy, but I'm having a hard time figuring out how to recode some data in R. I have a variable "numpeers" which is valued 1, 2, or 3. I also have three other variables called "w9zd9_1," "w9zd9_2," and "w9zd9_3." I want to use these variables to create a new item called "distot." Specifically, here is what I want to do: -if numpeers=1, then distot=w9zd9_1 -if numpeers=2, then distot=(w9zd9_1 + w9zd9_2)/2 -if numpeers=3, then distot=(w9zd9_1 + w9zd9_2 + w9zd9_3)/3 Thanks, in advance, for the help! Best, -c -- Casey A. Klofstad University of Miami Department of Political Science Coral Gables, FL klofstad(a)gmail.com http://www.as.miami.edu/personal/cklofstad/ - Zelig Mailing List, served by Harvard-MIT Data Center Send messages: zelig(a)lists.gking.harvard.edu [un]subscribe Options: http://lists.gking.harvard.edu/?info=zelig Zelig program information: http://gking.harvard.edu/zelig/

14 years, 7 months

Re: [matchit] Re: [zelig] full matching

by Peter Steiner

Hi Liz, thanks for your references. I tried to use the matchit() weights, but the weights either didn't change results very much or sim() produced an error message (as described below). Thanks, Peter Stuart, Elizabeth A. wrote: > Hi Peter, > Kosuke is right; full matching essentially uses all observations and > so your outcome models have to take into account the structure of the > matching in some way, either by using the subclasses that are formed > (e.g., through random or fixed effects) or by using the weights that > matchit produces (which are essentially frequency weights). I am not > sure your zelig code below does either of those things; it looks like > it is just taking a simple difference in means of the outcome in the > matched samples (which will be the same as the original samples). > > Ben Hansen discusses the subclass-specific fixed effects approach in > his 2004 paper (http://www.stat.lsa.umich.edu/%7Ebbh/hansen2004.pdf), > <http://www.stat.lsa.umich.edu/%7Ebbh/hansen2004.pdf%29,> or you can > also see description of both approaches in my paper on full matching > (Stuart, E.A., and Green, K.M. (2008). Using Full Matching to Estimate > Causal Effects in Non-Experimental Studies: Examining the Relationship > between Adolescent Marijuana Use and Adult Outcomes. /Developmental > Psychology/ 44(2): 395-406.) > > All of that said, I do also agree with Kosuke that more work is needed > to better understand the best approaches for doing analyses following > full matching. But I hope this helps. > > Thanks, > Liz > > On 9/11/09 9:55 AM, "Kosuke Imai" <kimai(a)Princeton.EDU> wrote: > > Full matching tries to use as many observations as possible, and > so if you > just do diff-in-means you will get a poor estimate (this is close > to not > doing any adjustment). I think that it is an open research > question as to > how to do model-based adjustment after full matching, and so I > don't want > to make any strong suggestion here. One way to do it however is > to use > either fixed or random effect type models so that the strata > created by > full matching can be incorporated into the model. But, others may > have > better ideas. > > Kosuke > > -- > Department of Politics > Princeton University > http://imai.princeton.edu > > On Wed, 9 Sep 2009, Peter Steiner wrote: > > > Hi, > > > > after briefly outlining my problem with the full-matching option > in Zelig and > > MatchIt to Imai, I post it in more detail and with the LaLonde data. > > Unfortunately, I was not able to figure out how to correctly > implement > > full-matching with Zelig using all your documentation. > > Moreover, the help in R seems to be sometimes inconsistent (the > help file for > > match.data() differs for the MatchIt and Zelig Package; one > suggests that > > weights are allowed (MatchIt) the other does not show an option > for weights > > (Zelig)) > > > > However, I'm interested in estimating ATE using full matching. A > > straightforward way seems to be the following: > > > > m.out <- matchit(treat ~ age + educ + black + hispan + > > nodegree + married + re74 + re75, data = lalonde, > > method = 'full', discard = 'both') > > z.out1 <- zelig(re78 ~ age + educ + black + hispan + nodegree + > > married + re74 + re75, data = match.data(m.out, "control"), > > model = 'ls') > > z.out0 <- zelig(re78 ~ age + educ + black + hispan + nodegree + > > married + re74 + re75, data = match.data(m.out, "treat"), > > model = 'ls') > > x.out1 <- setx(z.out1, data = match.data(m.out, "treat"), cond = T) > > x.out0 <- setx(z.out0, data = match.data(m.out, "control"), cond = T) > > s.out1 <- sim(z.out1, x = x.out1) > > s.out0 <- sim(z.out0, x = x.out0) > > ate.all <- c(s.out1$qi$att.ev, -s.out0$qi$att.ev) > > c(mean(ate.all), sd(ate.all)) > > > > This results in an ATE of 1335.5, which seems at least plausible. > But here is > > my problem. Using only full-matching without additional > covariance adjustment > > by setting > > > > z.out1 <- zelig(re78 ~ 1, data = match.data(m.out, "control"), > > model = 'ls') > > z.out0 <- zelig(re78 ~ 1, data = match.data(m.out, "treat"), > > model = 'ls') > > > > results in an ATE of 86.88, indicating a very poor performance of > full > > matching. (Using Ben Hansen's R package alone, full matching > would reduce > > almost all bias.) So I tried to use the weights generated by the > matchit() > > command. Using weights in match.data(..., weights = 'weights') > doesn't change > > the estimates very much. And whenever I include weights in > zelig() the sim() > > command produces an Error. > > > > So my question is, what is the correct way to perform > full-matching in Zelig? > > > > I have another question, concerning nearest neighbor matching > (1:1 matching, > > for instance). What does ATE exactly estimate in Zelig? ATE of > the original > > (unmatched) sample or ATE of the matched sample only (note that > the covariate > > distribution of the original and matched sample usually differ > considerably). > > I assume the latter. > > > > Thanks a lot, > > Peter > > > > > > -- > MatchIt mailing list served by Harvard-MIT Data Center > List Address: matchit(a)lists.gking.harvard.edu > Subscribe/Unsubscribe: http://lists.gking.harvard.edu/?info=matchit > MatchIt Software and Documentation: http://gking.harvard.edu/matchit/ > Browse/Search List Archive: > http://lists.hmdc.harvard.edu/lists/matchit/ > > -- : Peter M. Steiner : Northwestern University : Institute for Policy Research : 2040 Sheridan Road : Evanston, IL 60208 : p ++1 847 467 1479 : p-steiner(a)northwestern.edu - Zelig Mailing List, served by Harvard-MIT Data Center Send messages: zelig(a)lists.gking.harvard.edu [un]subscribe Options: http://lists.gking.harvard.edu/?info=zelig Zelig program information: http://gking.harvard.edu/zelig/

14 years, 7 months

full matching

by Peter Steiner

Hi, after briefly outlining my problem with the full-matching option in Zelig and MatchIt to Imai, I post it in more detail and with the LaLonde data. Unfortunately, I was not able to figure out how to correctly implement full-matching with Zelig using all your documentation. Moreover, the help in R seems to be sometimes inconsistent (the help file for match.data() differs for the MatchIt and Zelig Package; one suggests that weights are allowed (MatchIt) the other does not show an option for weights (Zelig)) However, I'm interested in estimating ATE using full matching. A straightforward way seems to be the following: m.out <- matchit(treat ~ age + educ + black + hispan + nodegree + married + re74 + re75, data = lalonde, method = 'full', discard = 'both') z.out1 <- zelig(re78 ~ age + educ + black + hispan + nodegree + married + re74 + re75, data = match.data(m.out, "control"), model = 'ls') z.out0 <- zelig(re78 ~ age + educ + black + hispan + nodegree + married + re74 + re75, data = match.data(m.out, "treat"), model = 'ls') x.out1 <- setx(z.out1, data = match.data(m.out, "treat"), cond = T) x.out0 <- setx(z.out0, data = match.data(m.out, "control"), cond = T) s.out1 <- sim(z.out1, x = x.out1) s.out0 <- sim(z.out0, x = x.out0) ate.all <- c(s.out1$qi$att.ev, -s.out0$qi$att.ev) c(mean(ate.all), sd(ate.all)) This results in an ATE of 1335.5, which seems at least plausible. But here is my problem. Using only full-matching without additional covariance adjustment by setting z.out1 <- zelig(re78 ~ 1, data = match.data(m.out, "control"), model = 'ls') z.out0 <- zelig(re78 ~ 1, data = match.data(m.out, "treat"), model = 'ls') results in an ATE of 86.88, indicating a very poor performance of full matching. (Using Ben Hansen's R package alone, full matching would reduce almost all bias.) So I tried to use the weights generated by the matchit() command. Using weights in match.data(..., weights = 'weights') doesn't change the estimates very much. And whenever I include weights in zelig() the sim() command produces an Error. So my question is, what is the correct way to perform full-matching in Zelig? I have another question, concerning nearest neighbor matching (1:1 matching, for instance). What does ATE exactly estimate in Zelig? ATE of the original (unmatched) sample or ATE of the matched sample only (note that the covariate distribution of the original and matched sample usually differ considerably). I assume the latter. Thanks a lot, Peter -- : Peter M. Steiner : Northwestern University : Institute for Policy Research : 2040 Sheridan Road : Evanston, IL 60208 : p ++1 847 467 1479 : p-steiner(a)northwestern.edu - Zelig Mailing List, served by Harvard-MIT Data Center Send messages: zelig(a)lists.gking.harvard.edu [un]subscribe Options: http://lists.gking.harvard.edu/?info=zelig Zelig program information: http://gking.harvard.edu/zelig/

14 years, 7 months

estimation uncertainty for random effects in mixed models?

by Malcolm Fairbrother

Dear all, This is my first post here--hope this is appropriate. I'm trying to figure out how to simulate quantities of interest from a mixed (multilevel) logit model. The logit.mixed documentation makes this look straightforward, but I'd like to know what's going on "under the hood". Accounting for estimation uncertainty for the fixed effects is straightforward--just the same as a single-level model. But what about taking into account uncertainty about the random effects? I've searched around on the Internet (and read the lme4 documentation) and haven't had any luck figuring out how to do this in R. (I realize that the random effects at each level are by definition a distribution with a mean of 0, so I guess in many instances simulations for quantities of interest can simply leave out the random effects, because their mean effect is 0. But if we want to simulate for a particular set of X values in a specific group versus in another group, the random effect will clearly matter.) Doug Bates is on record as saying that standard errors don't make sense for the random effects (http://www.stat.columbia.edu/~cook/movabletype/archives/2008/12/uncertainty… ), since each estimate is not likely to follow a symmetrical distribution. And lme4's mcmcsamp(), as I understand it, is not operational. Gelman and Hill (2007) suggest using a Bayesian approach and BUGS--getting a posterior distribution of estimates--but I don't know BUGS and would rather not have to learn just for this. MLwiN includes SEs for random effects, but doesn't have a straightforward facility for estimating quantities of interest. So my question (if this all makes sense) is how Zelig deals with this issue. Is it finding some way to take into account estimation uncertainty for the random effects? Or just treating each random effect as a given, and only taking into account estimation uncertainty for the fixed effects? Any illumination about this point would be much appreciated. Thanks, Malcolm Fairbrother University of Bristol

14 years, 7 months

Re: [zelig] zelig - Error in inherits(x, "factor")

by Olivia Lau

Hi, Zebedee, It seems that you're missing some important information in the zelig call. The most important thing for the zelig function is that you have all the variables in one data frame, which is specified using the data argument in the zelig call. To insert variables in your data frame Data, the line dtrdseas.log.CHPK.ts <- log(CHPK.ts)-(lm(log(CHPK.ts)~t+t2+t3+t3+sin.t+cos.t)$fit) should be Data$dtrdseas.log.CHPK.ts<- log(CHPK.ts) - (lm(log(CHPK.ts)~t+t2+t3+t3+sin.t+cos.t)$fit) where the $ inserts this variable in the Data dataframe. Next, the zelig call should be: z.out <- zelig(Diff(dtrdseas.log.CHPK.ts, 0) ~ lag.y(3) + lag.eps(3), model = "arima", data = Data) where Data is the data frame into which you have inserted your dtrdseas.log.CHPK.ts variable. Please let us know whether this fixes the problem. Yours, Olivia Lau Zelig[[3]] On Sat, Sep 5, 2009 at 10:37 AM, Zebedee Nii-Naate<zebedeeniinaate(a)gmail.com> wrote: > Thank you Olivia. > > I'm using R version 2.9.2 and Zelig version 3.4-5 (built: 2009-03-13). I'm > praying your colleague is able to have a look early week, as I'm due to hand > in my dissertation on Wed of next week. > > btw - I think you forgot to cc in colleague. > > Kind regards, > > Zebedee > > 2009/9/5 Olivia Lau <olivia.lau(a)post.harvard.edu> >> >> Dear Zebedee, >> >> The student who worked on ARIMA is currently at APSA. I've cc'ed him >> above, and hopefully he'll have a chance to take a look when he gets >> back next week. >> >> Which version of Zelig and which version of R are you using? >> >> Yours, >> Olivia >> >> On Sat, Sep 5, 2009 at 3:42 AM, Zebedee >> Nii-Naate<zebedeeniinaate(a)gmail.com> wrote: >> > Dear all, >> > >> > >> > >> > Good morning all. I’m trying to us the Zelig::zelig function to >> > estimate an >> > ARIMA model after I have removed the seasonal and trend components of a >> > time >> > series data (CHPH.ts) that I have imported from a CV file. Please see >> > syntax below. >> > >> > >> > >> > All seems to be ok apart from when I run: z.out <- >> > zelig(Diff(dtrdseas.log.CHPK.ts, 0) ~ lag.y(3) + lag.eps(3), model = >> > "arima" >> > , ,) >> > >> > >> > >> > I get the following error message: >> > >> > >> > >> > Error in inherits(x, "factor") : element 1 is empty; >> > >> > the part of the args list of 'class' being evaluated was: >> > >> > (data) >> > >> > >> > >> > I’ve tried changing the syntax to include: data = Data and data = >> > dtrdseas.log.CHPK.ts plus many alternatives. Nothing seems to work. >> > I’d be >> > very grateful for advice on the error message and how force the zelig >> > function to read dtrdseas.log.CHPK.ts in the zelig function. >> > >> > >> > >> > I’ve read material on the zelig website (http://gking.harvard.edu/zelig) >> > but >> > I’m still lost. >> > >> > >> > >> > Syntax >> > >> > #Import Data from csv file >> > >> > Data >> > >> > <-read.csv("C:/Users/Zebedee/Documents/Dissertation/Data/Germany.csv",header=T) >> > >> > >> > >> > #Name the first column of csv file as Data and format the Data >> > >> > Datecol<-Data[,1] >> > >> > as.Date(as.character(Datecol),format="%Y%m%d") >> > >> > >> > >> > #Report the dimensions of Data and list >> > >> > dim(Data) >> > >> > names (Data) >> > >> > >> > >> > ## DESCRIPTIVE STATISTICS >> > >> > # Plots in levels >> > >> > CH.ts <- ts(Data[,2], start = c(2002, 1), end = c(2008, 22), freq = 52) >> > >> > PK.ts <- ts(Data[,3], start = c(2002, 1), end = c(2008, 22), freq = 52) >> > >> > BF.ts <- ts(Data[,4], start = c(2002, 1), end = c(2008, 22), freq = 52) >> > >> > CHPK.ts<- ts(Data[,5], start = c(2002, 1), end = c(2008, 22), freq = 52) >> > >> > CHBF.ts<- ts(Data[,6], start = c(2002, 1), end = c(2008, 22), freq = 52) >> > >> > >> > >> > #Deseasonalise and detrend data >> > >> > t<-1: length(log(CHPK.ts)) >> > >> > t2<-t^2 >> > >> > t3<-t^3 >> > >> > sin.t<-sin(2*pi*t/52) >> > >> > cos.t<-cos(2*pi*t/52) >> > >> > >> > >> > par(mfrow=c(2,1)) >> > >> > plot(log(CHPK.ts)) >> > >> > trd.log.CHPK.ts<-lm(log(CHPK.ts)~t+t2+t3+sin.t+cos.t) >> > >> > lines(lm(log(CHPK.ts)~t+t2+t3+sin.t+cos.t)$fit,col=2,lwd=2) >> > >> > summary(lm(log(CHPK.ts)~t+t2+t3+sin.t+cos.t)) >> > >> > AIC(trd.log.CHPK.ts) >> > >> > tsDyn::BIC(trd.log.CHPK.ts) >> > >> > lmtest::bptest(trd.log.CHPK.ts) >> > >> > lmtest::resettest(trd.log.CHPK.ts) >> > >> > lmtest::bgtest(trd.log.CHPK.ts, order = 52) >> > >> > lmtest::dwtest(trd.log.CHPK.ts) >> > >> > dtrdseas.log.CHPK.ts<- log(CHPK.ts) - >> > (lm(log(CHPK.ts)~t+t2+t3+t3+sin.t+cos.t)$fit) >> > >> > plot(dtrdseas.log.CHPK.ts) >> > >> > >> > >> > #Estimate ARIMA model >> > >> > library(Zelig) >> > >> > z.out <- zelig(Diff(dtrdseas.log.CHPK.ts, 0) ~ lag.y(3) + lag.eps(3), >> > model >> > = "arima" , ,) >> > >> > summary(z.out) >> > >> > x.out <- setx(z.out) >> > >> > s.out <- sim(z.out) >> > >> > s.out <- sim(z.out, x = x.out) >> > >> > >> > >> > Kind regards, >> > >> > >> > >> > Zebedee >> > >> > > > - Zelig Mailing List, served by Harvard-MIT Data Center Send messages: zelig(a)lists.gking.harvard.edu [un]subscribe Options: http://lists.gking.harvard.edu/?info=zelig Zelig program information: http://gking.harvard.edu/zelig/

14 years, 8 months

zelig - Error in inherits(x, "factor")

by Zebedee Nii-Naate

Dear all, Good morning all. I'm trying to us the Zelig::zelig function to estimate an ARIMA model after I have removed the seasonal and trend components of a time series data (CHPH.ts) that I have imported from a CV file. Please see syntax below. All seems to be ok apart from when I run: z.out <- zelig(Diff(dtrdseas.log.CHPK.ts, 0) ~ lag.y(3) + lag.eps(3), model = "arima" , ,) I get the following error message: Error in inherits(x, "factor") : element 1 is empty; the part of the args list of 'class' being evaluated was: (data) I've tried changing the syntax to include: data = Data and data = dtrdseas.log.CHPK.ts plus many alternatives. Nothing seems to work. I'd be very grateful for advice on the error message and how force the zelig function to read dtrdseas.log.CHPK.ts in the zelig function. I've read material on the zelig website (http://gking.harvard.edu/zelig) but I'm still lost Syntax #Import Data from csv file Data <-read.csv("C:/Users/Zebedee/Documents/Dissertation/Data/Germany.csv",header =T) #Name the first column of csv file as Data and format the Data Datecol<-Data[,1] as.Date(as.character(Datecol),format="%Y%m%d") #Report the dimensions of Data and list dim(Data) names (Data) ## DESCRIPTIVE STATISTICS # Plots in levels CH.ts <- ts(Data[,2], start = c(2002, 1), end = c(2008, 22), freq = 52) PK.ts <- ts(Data[,3], start = c(2002, 1), end = c(2008, 22), freq = 52) BF.ts <- ts(Data[,4], start = c(2002, 1), end = c(2008, 22), freq = 52) CHPK.ts<- ts(Data[,5], start = c(2002, 1), end = c(2008, 22), freq = 52) CHBF.ts<- ts(Data[,6], start = c(2002, 1), end = c(2008, 22), freq = 52) #Deseasonalise and detrend data t<-1: length(log(CHPK.ts)) t2<-t^2 t3<-t^3 sin.t<-sin(2*pi*t/52) cos.t<-cos(2*pi*t/52) par(mfrow=c(2,1)) plot(log(CHPK.ts)) trd.log.CHPK.ts<-lm(log(CHPK.ts)~t+t2+t3+sin.t+cos.t) lines(lm(log(CHPK.ts)~t+t2+t3+sin.t+cos.t)$fit,col=2,lwd=2) summary(lm(log(CHPK.ts)~t+t2+t3+sin.t+cos.t)) AIC(trd.log.CHPK.ts) tsDyn::BIC(trd.log.CHPK.ts) lmtest::bptest(trd.log.CHPK.ts) lmtest::resettest(trd.log.CHPK.ts) lmtest::bgtest(trd.log.CHPK.ts, order = 52) lmtest::dwtest(trd.log.CHPK.ts) dtrdseas.log.CHPK.ts<- log(CHPK.ts) - (lm(log(CHPK.ts)~t+t2+t3+t3+sin.t+cos.t)$fit) plot(dtrdseas.log.CHPK.ts) #Estimate ARIMA model library(Zelig) z.out <- zelig(Diff(dtrdseas.log.CHPK.ts, 0) ~ lag.y(3) + lag.eps(3), model = "arima" , ,) summary(z.out) x.out <- setx(z.out) s.out <- sim(z.out) s.out <- sim(z.out, x = x.out) Kind regards, Zebedee

14 years, 8 months

zelig - Error in inherits(x, "factor")

by Zebedee Nii-Naate

Dear all, Good morning all. I'm trying to us the Zelig::zelig function to estimate an ARIMA model after I have removed the seasonal and trend components of a time series data (CHPH.ts) that I have imported from a CV file. Please see syntax below. All seems to be ok apart from when I run: z.out <- zelig(Diff(dtrdseas.log.CHPK.ts, 0) ~ lag.y(3) + lag.eps(3), model = "arima" , ,) I get the following error message: Error in inherits(x, "factor") : element 1 is empty; the part of the args list of 'class' being evaluated was: (data) I've tried changing the syntax to include: data = Data and data = dtrdseas.log.CHPK.ts plus many alternatives. Nothing seems to work. I'd be very grateful for advice on the error message and how force the zelig function to read dtrdseas.log.CHPK.ts in the zelig function. I've read material on the zelig website (http://gking.harvard.edu/zelig) but I'm still lost. Syntax #Import Data from csv file Data <-read.csv("C:/Users/Zebedee/Documents/Dissertation/Data/Germany.csv",header =T) #Name the first column of csv file as Data and format the Data Datecol<-Data[,1] as.Date(as.character(Datecol),format="%Y%m%d") #Report the dimensions of Data and list dim(Data) names (Data) ## DESCRIPTIVE STATISTICS # Plots in levels CH.ts <- ts(Data[,2], start = c(2002, 1), end = c(2008, 22), freq = 52) PK.ts <- ts(Data[,3], start = c(2002, 1), end = c(2008, 22), freq = 52) BF.ts <- ts(Data[,4], start = c(2002, 1), end = c(2008, 22), freq = 52) CHPK.ts<- ts(Data[,5], start = c(2002, 1), end = c(2008, 22), freq = 52) CHBF.ts<- ts(Data[,6], start = c(2002, 1), end = c(2008, 22), freq = 52) #Deseasonalise and detrend data t<-1: length(log(CHPK.ts)) t2<-t^2 t3<-t^3 sin.t<-sin(2*pi*t/52) cos.t<-cos(2*pi*t/52) par(mfrow=c(2,1)) plot(log(CHPK.ts)) trd.log.CHPK.ts<-lm(log(CHPK.ts)~t+t2+t3+sin.t+cos.t) lines(lm(log(CHPK.ts)~t+t2+t3+sin.t+cos.t)$fit,col=2,lwd=2) summary(lm(log(CHPK.ts)~t+t2+t3+sin.t+cos.t)) AIC(trd.log.CHPK.ts) tsDyn::BIC(trd.log.CHPK.ts) lmtest::bptest(trd.log.CHPK.ts) lmtest::resettest(trd.log.CHPK.ts) lmtest::bgtest(trd.log.CHPK.ts, order = 52) lmtest::dwtest(trd.log.CHPK.ts) dtrdseas.log.CHPK.ts<- log(CHPK.ts) - (lm(log(CHPK.ts)~t+t2+t3+t3+sin.t+cos.t)$fit) plot(dtrdseas.log.CHPK.ts) #Estimate ARIMA model library(Zelig) z.out <- zelig(Diff(dtrdseas.log.CHPK.ts, 0) ~ lag.y(3) + lag.eps(3), model = "arima" , ,) summary(z.out) x.out <- setx(z.out) s.out <- sim(z.out) s.out <- sim(z.out, x = x.out) Kind regards, Zebedee

14 years, 8 months

2024

2023

2022

2021

2020

2019

2018

2017

2016

2015

2014

2013

2012

2011

2010

2009

2008

2007

2006

2005

2004

2003

Zelig September 2009