Cointegration in R

Essentially, it seeks to find stationary linear combinations of the two vectors
test two series for integration and return the p-value indicating the likelihood of correlation

This runs an augmented Dickey-Fuller test and will return a p-value indicating whether the series are mean-reverting or not. You can use the typical p-value as a test of significance if you like(ie, a p-value below .05 indicates a mean-reverting spread), or you can use an alternate value. This assumes that your two series were observed at the same time points.

gld <- read.csv("http://ichart.finance.yahoo.com/table.csv?s=6947.KL&ignore=.csv", stringsAsFactors=F)
gdx <- read.csv("http://ichart.finance.yahoo.com/table.csv?s=6012.KL&ignore=.csv", stringsAsFactors=F)

Take the intersection of the two zoo objects.  That will create one zoo object with the observations common to both datasets.

#The seventh column contains the adjusted close.
#The first column contains dates.

gld <- zoo(gld[,7], as.Date(gld[,1]))
gdx <- zoo(gdx[,7], as.Date(gdx[,1]))

The merge function can combine two zoo objects. so we merge them.
t.zoo <- merge(gld, gdx, all=FALSE)

# At this point, t.zoo is a zoo object with two columns: gld and gdx.
# Most statistical functions expect a data frame for input,
# so we create a data frame here.
#
t <- as.data.frame(t.zoo)


First we construct the spread, then we test the spread for a unit root. 

It the spread has a root inside the unit circle, the underlying securities are cointegrated.

The lm function builds linear regression models using  ordinary least squares(OLS).
# We build the linear model, m, forcing a zero intercept,
# then we extract the model's first regression coefficient.
m <- lm(gld ~ gdx + 0, data=t)
beta <- coef(m)[1]
cat("Assumed hedge ratio is", beta, "\n")

sprd <- t$gld - beta*t$gdx

The Augmented Dickey-Fuller test is a basic statistical test for a unit root, and several R packages implement that test.  Here, we will use the adf.test function which is implemented in the tseries package.  The function returns an object which contains the test results. In particular, it contains the p-value that we want.
library(tseries)            # Load the tseries package

# Setting alternative="stationary" chooses the appropriate test.
# Setting k=0 forces a basic (not augmented) test. 
ht <- adf.test(sprd, alternative="stationary", k=0)
cat("ADF p-value is", ht$p-value, "\n")

The adf.test function essentially detrends your data before testing for stationarity.
If your data contains a strong trend consider the  fUnitRoots package which contains the adfTest function

if (ht$p.value < 0.05) {
    cat("The spread is likely mean-reverting.\n")
} else {
    cat("The spread is not mean-reverting.\n")
}

gld <- read.csv("http://ichart.finance.yahoo.com/table.csv?s=6947.KL&ignore=.csv", stringsAsFactors=F)

SUMMARY


gld <- read.csv("http://ichart.finance.yahoo.com/table.csv?s=6012.KL&ignore=.csv", stringsAsFactors=F)
gdx <- read.csv("http://ichart.finance.yahoo.com/table.csv?s=4863.KL&ignore=.csv", stringsAsFactors=F)
gld <- zoo(gld[,7], as.Date(gld[,1]))
gdx <- zoo(gdx[,7], as.Date(gdx[,1]))
t.zoo <- merge(gld, gdx, all=FALSE)
t <- as.data.frame(t.zoo)
m <- lm(gld ~ gdx + 0, data=t)
beta <- coef(m)[1]
sprd <- t$gld - beta*t$gdx
ht <- adf.test(sprd, alternative="stationary", k=0)

if (ht$p.value < 0.05) {
     cat("The spread is likely mean-reverting.\n")
 } else {
     cat("The spread is not mean-reverting.\n")
 }(sprd, alternative="stationary", k=0)

Singular spectrum analysis basics in R

Using the following R package
 http://cran.r-project.org/web/packages/Rssa/index.html

and my time series as the variable: z.ts

Summary: Singular spectrum analysis for time series
Anatoly Zhigljavsky

Singular Spectrum Analysis a technique of times series analysis and forecasting.
Aim is to decompose the original series into a sum of smaller number of interpretable components such as:
slowly varying trend, oscillatory components and a structureless noise.

based on singular value decomposition (SVD) of a specific matrix constructed upon time series.

SSA is a model-free technique because no assumptions such as parametric model nor stationary type

condition is required.

Basic SSA
X= construct the trajectory matrix (lagged vectors)
this matrix is a Hankel Matrix, all elements along the diagonal are equal

the SVD of matrix XX(transpose) yields a collection of L eigenvalues and eigenVectors.

Basic SSA can be used for smoothing,filtration, noise reduction, extraction of trends of different

resolution, extraction of periodicities in the form of modulated harmonics, gap-filling

One of the requirements of SSA is a continuous time series with no Gaps.





in R statistics
s <- new.ssa(z.ts)

suitable grouping of the elementary time series is required via looking at the eigenplots of the

decomposition

plot(s, type = "series", groups = list(1:4)) #Plot the first 4 reconstructed components
plot(s, type = "values") #Plot the eigenvalues

examine the so-called w-correlation matrix
# Calculate the w-correlation matrix between first 10 series
w <- wcor(s, groups = 1:10)
print(w)
plot(w)

reconstruction of the time-series using the selected grouping
# Reconstruct the series, grouping elementary series 2, 3 and 4, 5.
r <- reconstruct(s, groups = list(1, c(2,3), c(4,5)))
plot(r$F1, col = "black")
lines(r$F1 + r$F2, col = "red")
lines(r$F1 + r$F2 + r$F3, col = "blue")

Summary of steps to get the graph

s <- new.ssa(z.ts) # Perform the decomposition using the default window length
summary(s)        # Show various information about the decomposition
plot(s)           # Show the plot of the eigenvalues
f <- reconstruct(s, groups = list(1, c(2, 3), 4)) # Reconstruct into 3 series
plot(z.ts)         # Plot the original series
lines(f$F1, col = "blue")            # Extract the trend
lines(f$F1+f$F2, col = "red")        # Add the periodicity
lines(f$F1+f$F2+f$F3, col = "green") # Add slow-varying component

Forcast
# Produce 5 forecasted values and confidence bounds of the series using
# the first 3 eigentriples as a base space for the forecast.
bforecast(s, group = 1:3, len = 5)

Using R as a Tool and getting your dataset

Download R
statistical computing and graphics tool.
www.r-project.org/

Download R Packages from http://cran.r-project.org/web/package/available_packages_by_name.html
forcast package with these dependencies
tseries, fracdiff, zoo, Rcpp (≥ 0.9.10), RcppArmadillo (≥ 0.2.35)

Download RStudio
http://rstudio.org/
An IDE for R use RStudio interface to download and install the packages (See image)

Find your Stock Symbol
http://finance.yahoo.com/
Search by entering the name and click "Get Quotes"

Download your Dataset using this URL (change to your desired stock symbol)
http://ichart.finance.yahoo.com/table.csv?s=4677.KL&a=0&b=1&c=2000&d=08&e=9&f=2012&g=d&ignore=.csv

s=4677.KL : my stock code
after s= goes the ticker symbol, after a= the start month (minus 1), after b= the start day, c= the start year and so on. The final g= parameter lets you choose between getting historical stock information on a daily, weekly, or monthly basis.

You have your Tool and Data. Lets Try it (Variables are CASE SENSITIVE):

Rename the download CSV file to the name of the stock, in my case table.csv->YTL.csv
Import the data to RStudio
YTL <- read.csv("C:/Users/user/Desktop/YTL.csv")

Convert the dataset to a timeseries
YTL <- ts(YTL[,-1], start=2000, frequency=4)

start is the year we are observing the data
frequncy is the number of observations per unit of time.

Plot and View the Summary
plot(YTL)
summary(YTL)

Any Correlation between Volume and Closing Price
cor.test(YTL[,"Volume"],YTL[,"Close"])

Scatterplot matrix of the variables:
 pairs(as.data.frame(YTL))

Monthly Plot:
monthplot(YTL[,"Close"])
monthplot(YTL[,"Volume"])