Posts Tagged ‘R’

Slides of Bass Diffusion Model

Published by chengjun on April 26th, 2012

http://weblab.com.cityu.edu.hk/blog/chengjun/files/2012/04/Bass-Diffusion-Model.pdf

Introduction

I had made slides to understand the bass diffusion model which is proposed by Frank M. Bass in 1969. The model proposes that diffusion is motivated both by inovativeness and imitation: first, the innovative early adopters adopted the products, after which the followers will imitate them and adopt the products.

•The Bass model coefficient (parameter) of innovation is p.
•The Bass model coefficient (parameter) of imitation is q.

The slides will show you how to deprive the equation of Bass diffusion model by both discrete-time model and continuous model (using Hazard rate).
“The probability of adopting by those who have not yet adopted is a linear function of those who had previously adopted.”
By solving the differential equation of bass diffusion model using Mathematica and Simulating the model using R code.
# basss diffusion model
# chengjun, 20120424@canberra

# refer to http://en.wikipedia.org/wiki/Bass_diffusion_model
# and http://book.douban.com/subject/4175572/discussion/45634092/

# BASS Diffusion Model
# three parameters:
# the total number of people who eventually buy the product, m;
# the coefficient of innovation, p;
# and the coefficient of imitation, q

# example
T79 <- 1:10
Tdelt <- (1:100) / 10
Sales <- c(840,1470,2110,4000, 7590, 10950, 10530, 9470, 7790, 5890)
Cusales <- cumsum(Sales)
Bass.nls <- nls(Sales ~ M * ( ((P+Q)^2 / P) * exp(-(P+Q) * T79) ) /(1+(Q/P)*exp(-(P+Q)*T79))^2,
                start = list(M=60630, P=0.03, Q=0.38))
summary(Bass.nls)

# get coefficient
Bcoef <- coef(Bass.nls)
m <- Bcoef[1]
p <- Bcoef[2]
q <- Bcoef[3]
# setting the starting value for M to the recorded total sales.
ngete <- exp(-(p+q) * Tdelt)

# plot pdf
Bpdf <- m * ( (p+q)^2 / p ) * ngete / (1 + (q/p) * ngete)^2
plot(Tdelt, Bpdf, xlab = "Year from 1979",ylab = "Sales per year", type='l')
points(T79, Sales)

# plot cdf
Bcdf <- m * (1 - ngete)/(1 + (q/p)*ngete)
plot(Tdelt, Bcdf, xlab = "Year from 1979",ylab = "Cumulative sales", type='l')
points(T79, Cusales)

# when q=0, only Innovator without immitators.
Ipdf<- m * ( (p+0)^2 / p ) * exp(-(p+0) * Tdelt) / (1 + (0/p) * exp(-(p+0) * Tdelt))^2
# plot(Tdelt, Ipdf, xlab = "Year from 1979",ylab = "Isales per year", type='l')
Impdf<-Bpdf-Ipdf
plot(Tdelt, Bpdf, xlab = "Year from 1979",ylab = "Sales per year", type='l', col="red")
lines(Tdelt,Impdf,col="green")
lines(Tdelt,Ipdf,col="blue")

# when q=0, only Innovator without immitators.
Icdf<-m * (1 - exp(-(p+0) * Tdelt))/(1 + (0/p)*exp(-(p+0) * Tdelt))
# plot(Tdelt, Icdf, xlab = "Year from 1979",ylab = "ICumulative sales", type='l')
Imcdf<-m * (1 - ngete)/(1 + (q/p)*ngete)-Icdf
plot(Tdelt, Imcdf, xlab = "Year from 1979",ylab = "Cumulative sales", type='l', col="red")
lines(Tdelt,Bcdf,col="green")
lines(Tdelt,Icdf,col="blue")

Estimating Threshold of Time Series Using R

Published by chengjun on March 20th, 2012

Please note that I am NOT an expert in time series analysis. Therefore, I am not the ideal person to answer the technical questions on this topic. Please consider (1) raising your question on stackoverflow, (2) sending emails to the developer of related R packages, (3) joining related email groups, etc.

The method of estimating Threshold of Time Series Data has been developed by R. This post shows how to use the method by adopting two packages. First, I would like to highlight Bruce Hansen’s work in this field.

Bruce E. Hansen


Programs — Threshold Models


“Inference when a nuisance parameter is not identified under the null hypothesis.” Econometrica, (1996). [Download].

“Inference in TAR models.” Studies in Nonlinear Dynamics and Econometrics, (1997). [Download].

“Threshold effects in non-dynamic panels: Estimation, testing and inference.” Journal of Econometrics, (1999). [Download].

“Testing for Linearity.” Journal of Economic Surveys, (1999). [Download].

“Sample splitting and threshold estimation.” Econometrica, (2000). [Download].

“Threshold Autoregression with a Unit Root.” Econometrica (2001), with Mehmet Caner. [Download].

“How responsive are private transfers to income? Evidence from a laissez-faire economy.”
with Donald Cox and Emmanuel Jimenez, Journal of Public Economics, (2004), 88, 2193-2219. [Download].

“Testing for two-regime threshold cointegration in vector error correction models,” with Byeongseon Seo, Journal of Econometrics (2002). [Download].

“Instrumental Variable Estimation of a Threshold Model”, with Mehmet Caner, Econometric Theory, (2004), 20, 813-843. [Download].

# Chengjun WANG
# @anu, 20120320
#~~~~~~~~~~~~~~~~urca~~~~~~~~~~~~~~~~~~~~~~~~~#
# The first step is Unit root and cointegration Analysis (urca)
# install.packages("urca")
# load the package
library(urca)
# see the data which is used as an example
data(denmark)
data(finland)
#~~~~~~~~~~~~~~~threshold model using tsDyn~~~~~~~~~~~~~~~~~~~~~~~#
# install.packages("tsDyn")
# load the package of tsDyn
library(tsDyn)

# models in this package
availableModels()
 [1] "linear"  "nnetTs"  "setar"   "lstar"   "star"    "aar"     "lineVar"
 [8] "VECM"    "TVAR"    "TVECM"

#fit an AAR model:
mod #Summary data informations:
summary(mod)
#Diagnostic plots:
plot(mod)

# STAR model fitting with automatic selection of the number
# of regimes based on LM tests.

mod.star mod.star

addRegime(mod.star)

# Estimate a multivariate Threshold VAR
?TVAR  # ask r to introduce about TVAR

data(zeroyld)

data
TVAR(data, lag=2, nthresh=2, thDelay=1, trim=0.1, mTh=1, plot=TRUE)
TVAR.LRtest(data, lag=2, mTh=1,thDelay=1:2, nboot=3, plot=FALSE, trim=0.1, test="1vs")

# The one threshold (two regimes) gives a value of 10.698 for the
# threshold and 1 for the delay. Conditional on this values, the search
# for a second threshold (three regimes) gives 8.129. Starting from this
# values, a full grid search finds the same values and confims the first
# step estimation.

##simulate VAR as in Enders 2004, p 268
B1var1ts.plot(var1, type="l", col=c(1,2))

B2varcovvar2ts.plot(var2, type="l", col=c(1,2))

##Simulation of a TVAR with 1 threshold
Bsim
#estimate the new serie
TVAR(sim, lag=1, dummyToBothRegimes=TRUE)

##Bootstrap a TVAR with two threshold (three regimes)
data(zeroyld)
serieTVAR.sim(data=serie,nthresh=2, type="boot",mTh=1, Thresh=c(7,9))

##Check the bootstrap
cbind(TVAR.sim(data=serie,nthresh=2, type="check",mTh=1, Thresh=c(7,9)),serie)

# Estimate a Threshold Vector Error Correction model (VECM)
# Hansen, B. and Seo, B. (2002), Testing for two-regime threshold cointegration in vector error-correction models, Journal of Econometrics, 110, pages 293 - 318
data(zeroyld)
data
##Estimate a TVECM (we use here minimal grid, it should be usually much bigger!)

tv
print(tv)
summary(tv)

#Obtain diverse infos:
AIC(tv)
BIC(tv)

res.tv
#export the equations as Latex:
toLatex(tv)

How to Get Standard Regression Coefficient Using R

Published by admin on March 18th, 2012

Last week, Lexing showed us a picture of how different programming languages look like. Very Interesting. However, let’s dive to the question raised in the title of this post: How to Get Standard Regression Coefficient Using R?( In addition to that, I want to practice using SyntaxHighliter here.)

# e.g. I set up one regression
reg<-lm(view/top$DAYS~feo+fev+ffvv+frf+frfsm+frfgs
        +frfgvs+frflv+frfy+frfys+fvfa+fvfmd+fvocp+ov)
summary(reg)
# model diagnosis
layout(matrix(1:4,2,2))
plot(reg)

# according to the formula between standard regression coefficient
                beta.x=b.x*sd.x/sd.y
# calcuate the standard regression coefficient one by one
# for the beta of feo
b.x<-coef(reg)[2]
sd.y<-sd(view/top$DAYS)
sd.x<-sd(feo)
beta.x<-b.x*sd.x/sd.y
# it's dull to do it in this way!

# output standard coefficient using QuantPsyc library
library(QuantPsyc)# install.packages("QuantPsyc")
as.data.frame(lm.beta(reg))

Apparently, R is not well designed in this aspect, although we could get what we want, but it’s not so efficient and convenient compared with other commercial software. Hope this will be improved in the future.


Visualization of Dynamic Changes Using googleVis, r, and dropbox

Published by admin on September 12th, 2011

GoogleVis could be used to visualize the dynamic change of social pattern. Here I will test some examples.

library(googleVis)

data(package="googleVis")
# Data sets in package ‘googleVis’:
# Andrew            Hurricane Andrew: googleVis example data set
# CityPopularity    CityPopularity: googleVis example  data set
# Exports           Exports: googleVis example data set
# Fruits            Fruits: googleVis example data set
# OpenClose         OpenClose: googleVis example data set
# Population        Population: googleVis example data  set
# Regions           Regions: googleVis example data set
# Stock             Stock: googleVis example data set

# Visualizing the hurricane of Andrew

data(Andrew)
plot(Andrew)

AndrewGeoMap <- gvisGeoMap(Andrew, locationvar='LatLong', numvar='Speed_kt',
hovervar='Category',
options=list(width=600,height=300,
region='US', dataMode='Markers'))
plot(AndrewGeoMap)
AndrewGeoMap$html$chart
setwd("d:/r")
cat(AndrewGeoMap$html$chart, file="AndrewGeoMap.html")
# then you can find it in the work directory.
# Since WordPress doesn’t allow embedded JavaScript. I  put it into dropbox public folder and i past the

link is given here.

Further, I visualize the dynamic change of cell phone ratio in the world.

# fixed telephone use

ftel<-read.csv("d:/r/Fixed Telephone.csv", header = T, sep = ",", quote = "\"'",  dec = ".")[1:7263,]
ftelr<-read.csv("d:/r/Fixed Telephone ratio.csv", header = T, sep = ",", quote = "\"'",    dec = ".")

ftelm<-merge(ftelr, ftel, by=c("Country.or.Area","Year"))

ftelm<-cbind(ftelm[,1:3], ftelm[,5])
names(ftelm)<-c("Country or Area","Year","Percentage","Users")
library("googleVis")
fixed_telephone_motion<- gvisMotionChart(ftelm, idvar="Country or Area", timevar="Year", options=list(width=1024, height=768))plot(fixed_telephone_motion)
setwd("d:/r")
cat(fixed_telephone_motion$html$chart, file="fixed_telephone_motion.html")

The geographical distribution of cell phone usge (per 100 people)

The motion chart of cell phone usage (per 100 people)

# International news coverage visualization

Finally, i visualize the international news of different countries during August of 2011.

You may not be able to see the links above, I made a web page to show the result, click here:

http://weblab.com.cityu.edu.hk/chengjun/