Convert R Code to MATLAB for Detrended Cross Correlation Analysis

Click For Summary
SUMMARY

The discussion revolves around converting R code for Detrended Cross-Correlation Analysis (DCCA) into MATLAB. The R code utilizes libraries such as tseries, fracdiff, and matrixStats to compute beta coefficients and perform regression analysis on two time series. Key functions include DCCA_beta_avg, DCCA_beta_SE, and DCCA_beta_s, which handle various statistical calculations. Participants emphasize the importance of understanding the R code and its libraries before attempting the conversion to MATLAB.

PREREQUISITES
  • Familiarity with R programming language and syntax
  • Understanding of statistical concepts such as regression analysis and beta coefficients
  • Knowledge of MATLAB programming environment
  • Experience with time series analysis techniques
NEXT STEPS
  • Study the R libraries tseries, fracdiff, and matrixStats to understand their functions and parameters
  • Learn about the MATLAB equivalents for statistical functions used in the R code
  • Research the implementation of Detrended Cross-Correlation Analysis in MATLAB
  • Practice converting simple R functions to MATLAB to build familiarity with both languages
USEFUL FOR

Data scientists, statisticians, and programmers interested in time series analysis and those transitioning from R to MATLAB for statistical computations.

Tilfani
Messages
11
Reaction score
0

Homework Statement


I would like to convert this R code into matlab. This code perform a calculation of beta coefficient based on Detrended cross correlation analysis.

Homework Equations


Code: 
require(tseries)

require(fracdiff)

require(matrixStats)DCCA_beta_avg<-function(y,x,smin,smax,step){

  XTX<-var(x)*(length(x)-1)

  betas<-rep(0,(smax-smin)/step+1)

  for(s in seq(smin,smax,by=step)){

    betas[(s-smin)/step+1]<-DCCA_beta_sides(y,x,s)

  }

  DCCA_beta<-mean(betas)

  DCCA_res<-(y-DCCA_beta*x)-mean(y-DCCA_beta*x)

  DCCA_sigma2<-sum(DCCA_res^2)/(length(DCCA_res)-2)

  DCCA_SE<-sqrt(DCCA_sigma2/XTX)

  DCCA_R2<-1-var(DCCA_res)/var(y)

  OLS_beta<-lm(y~x)$coefficients[2]

  OLS_res<-(y-OLS_beta*x)-mean(y-OLS_beta*x)

  OLS_sigma2<-sum(OLS_res^2)/(length(OLS_res)-2)

  OLS_SE<-sqrt(OLS_sigma2/XTX)

  OLS_R2<-1-var(OLS_res)/var(y)

  return(c(OLS_beta,OLS_SE,OLS_R2,DCCA_beta,DCCA_SE,DCCA_R2))

}DCCA_beta<-DCCdccafunction(y,x,s){

  xx<-cumsum(x-mean(x))

  yy<-cumsum(y-mean(y))

  t<-1:length(xx)

  F2sj_xy<-runif(floor(length(xx)/s))

  F2sj_xx<-F2sj_xy

  for(ss in seq(1,(floor(length(xx)/s)*s),by=s)){

    F2sj_xy[(ss-1)/s+1]<-sum((summary(lm(xx[ss:(ss+s-1)]~t[ss:(ss+s-1)]))$residuals)*(summary(lm(yy[ss:(ss+s-1)]~t[ss:(ss+s-1)]))$residuals))/(s-1)

    F2sj_xx[(ss-1)/s+1]<-sum((summary(lm(xx[ss:(ss+s-1)]~t[ss:(ss+s-1)]))$residuals)*(summary(lm(xx[ss:(ss+s-1)]~t[ss:(ss+s-1)]))$residuals))/(s-1)

  }

  beta<-mean(F2sj_xy)/mean(F2sj_xx)

  return(beta)

}DCCA_beta_F<-function(y,x,s){

  xx<-cumsum(x-mean(x))

  yy<-cumsum(y-mean(y))

  t<-1:length(xx)

  F2sj_xy<-runif(floor(length(xx)/s))

  F2sj_xx<-F2sj_xy

  F2sj_yy<-F2sj_xy

  for(ss in seq(1,(floor(length(xx)/s)*s),by=s)){

    F2sj_xy[(ss-1)/s+1]<-sum((summary(lm(xx[ss:(ss+s-1)]~t[ss:(ss+s-1)]))$residuals)*(summary(lm(yy[ss:(ss+s-1)]~t[ss:(ss+s-1)]))$residuals))/(s-1)

    F2sj_xx[(ss-1)/s+1]<-sum((summary(lm(xx[ss:(ss+s-1)]~t[ss:(ss+s-1)]))$residuals)*(summary(lm(xx[ss:(ss+s-1)]~t[ss:(ss+s-1)]))$residuals))/(s-1)

    F2sj_yy[(ss-1)/s+1]<-sum((summary(lm(yy[ss:(ss+s-1)]~t[ss:(ss+s-1)]))$residuals)*(summary(lm(yy[ss:(ss+s-1)]~t[ss:(ss+s-1)]))$residuals))/(s-1)

  }

  beta<-mean(F2sj_xy)/mean(F2sj_xx)

  return(c(beta,mean(F2sj_xx),mean(F2sj_yy)))

  #return(c(beta,sum(F2sj_xx),sum(F2sj_yy)))

}DCCA_beta_SE<-function(y,x,s){

  r<-DCCA_beta_F(y,x,s)

  beta<-r[1]

  yhat<-beta*x

  alpha<-mean(y)-beta*mean(x)

  res<-y-yhat

  residuals<-res-mean(res)

  resres<-cumsum(residuals-mean(residuals))

  F2sj_res<-runif(floor(length(residuals)/s))

  t<-1:length(resres)

  for(ss in seq(1,(floor(length(residuals)/s)*s),by=s)){

    F2sj_res[(ss-1)/s+1]<-sum((summary(lm(resres[ss:(ss+s-1)]~t[ss:(ss+s-1)]))$residuals)*(summary(lm(resres[ss:(ss+s-1)]~t[ss:(ss+s-1)]))$residuals))/(s-1)

  }

  #SE<-mean(residuals^2)/((length(residuals)-2)*r[2])

  SE<-mean(F2sj_res)/((length(residuals)-2)*r[2])

  SE_a<-(mean(F2sj_res)/r[2])*(sum(x^2)/(length(residuals)*(length(residuals)-2)))

  R<-1-mean(F2sj_res)/(r[3])

  return(c(alpha,sqrt(SE_a),beta,sqrt(SE),R))

}DCCA_beta_SE_F<-function(y,x,s){

  r<-DCCA_beta_F(y,x,s)

  beta<-r[1]

  yhat<-beta*x

  alpha<-mean(y)-beta*mean(x)

  res<-y-yhat

  residuals<-res-mean(res)

  res_R<-y-x

  resres<-cumsum(residuals-mean(residuals))

  resres_R<-cumsum(res_R)

  F2sj_res<-runif(floor(length(residuals)/s))

  F2sj_res_R<-runif(floor(length(res_R)/s))

  t<-1:length(resres)

  for(ss in seq(1,(floor(length(residuals)/s)*s),by=s)){

    F2sj_res[(ss-1)/s+1]<-sum((summary(lm(resres[ss:(ss+s-1)]~t[ss:(ss+s-1)]))$residuals)*(summary(lm(resres[ss:(ss+s-1)]~t[ss:(ss+s-1)]))$residuals))/(s-1)

    F2sj_res_R[(ss-1)/s+1]<-sum((summary(lm(resres_R[ss:(ss+s-1)]~t[ss:(ss+s-1)]))$residuals)*(summary(lm(resres_R[ss:(ss+s-1)]~t[ss:(ss+s-1)]))$residuals))/(s-1)

  }

  #SE<-mean(residuals^2)/((length(residuals)-2)*r[2])

  #SE<-mean(F2sj_res)/((length(residuals)-2)*r[2])

  #SE<-mean(F2sj_res)/((length(F2sj_res)-2)*r[2]) #controlling for uncertainty connected to scales (higher scales have higher uncertainty due to lower number of blocks)

  SE<-mean(F2sj_res)/(ceiling(length(residuals)/s)*r[2]) #loosing d.f. due to fitting a and b in each box

  #SE_a<-(mean(F2sj_res)/r[2])*(sum(x^2)/(length(residuals)*(length(residuals)-2)))

  #SE_a<-(mean(F2sj_res)/r[2])*(sum(x^2)/(length(F2sj_res)*(length(F2sj_res)-2))) #controlling for uncertainty connected to scales (higher scales have higher uncertainty due to lower number of blocks)

  SE_a<-(mean(F2sj_res)/r[2])*(sum(x^2)/(length(residuals)*ceiling(length(residuals)/s))) #loosing d.f. due to fitting a and b in each box

  R<-1-mean(F2sj_res)/(r[3])

  #SSR_U<-sum(residuals^2)

  SSR_U<-sum(F2sj_res)

  #SSR_R<-sum((y-x)^2) #specific null: alpha=0, beta=1

  SSR_R<-sum(F2sj_res_R)

  #F_stat<-((SSR_R-SSR_U)/(SSR_U))*((length(residuals)-2)/2)

  #F_stat<-((SSR_R-SSR_U)/(SSR_U))*((length(F2sj_res)-2)/2) #controlling for uncertainty connected to scales (higher scales have higher uncertainty due to lower number of blocks)

  F_stat<-((SSR_R-SSR_U)/(SSR_U))*(ceiling(length(residuals)/s)/2) #loosing d.f. due to fitting a and b in each box

  F_p<-pf(F_stat,2,length(F2sj_res)-2,lower.tail=FALSE)

  return(c(alpha,sqrt(SE_a),beta,sqrt(SE),R,F_stat,F_p))

}DCCA_beta_s<-function(y,x,smin,smax,step){

  results<-matrix(rep(0,6*((smax-smin)/step+1)),ncol=6)

  for(s in seq(smin,smax,by=step)){

    beta<-DCCA_beta_SE(y,x,s)

    results[((s-smin)/step+1),1]<-s

    results[((s-smin)/step+1),2]<-beta[1]

    results[((s-smin)/step+1),3]<-beta[2]

    results[((s-smin)/step+1),4]<-beta[3]

    results[((s-smin)/step+1),5]<-beta[4]

    results[((s-smin)/step+1),6]<-beta[5]

  }

  return(results)

}DCCA_beta_s_F<-function(y,x,smin,smax,step){

  results<-matrix(rep(0,10*((smax-smin)/step+2)),ncol=10)

  for(s in seq(smin,smax,by=step)){

    beta<-DCCA_beta_SE_F(y,x,s)

    results[((s-smin)/step+1),1]<-s

    results[((s-smin)/step+1),2]<-beta[1]

    results[((s-smin)/step+1),3]<-beta[2]

    results[((s-smin)/step+1),4]<-2*pnorm(abs(beta[1]/beta[2]),lower.tail=FALSE)#p-value for null=0

    results[((s-smin)/step+1),5]<-beta[3]

    results[((s-smin)/step+1),6]<-beta[4]

    results[((s-smin)/step+1),7]<-2*pnorm(abs((beta[3]-1)/beta[4]),lower.tail=FALSE)#p-value for null=1

    results[((s-smin)/step+1),8]<-beta[5]

    results[((s-smin)/step+1),9]<-beta[6]

    results[((s-smin)/step+1),10]<-beta[7]

  }

  #results[(smax-smin)/step+2,2]<-mean(results[1:(dim(results)[1]-1),2])#A

  #results[(smax-smin)/step+2,5]<-mean(results[1:(dim(results)[1]-1),5])#B

  results[(smax-smin)/step+2,2]<-sum(results[1:(dim(results)[1]-1),2]*results[1:(dim(results)[1]-1),8])/sum(results[1:(dim(results)[1]-1),8])#A as R2(s) weighted

  results[(smax-smin)/step+2,5]<-sum(results[1:(dim(results)[1]-1),5]*results[1:(dim(results)[1]-1),8])/sum(results[1:(dim(results)[1]-1),8])#B as R2(s) weighted

  results[(smax-smin)/step+2,3]<-sqrt((sum(x^2)/length(x))*sum((y-results[(smax-smin)/step+2,2]-results[(smax-smin)/step+2,5]*x)^2)/((length(y)-dim(results)[1]+1)*sum((x-mean(x))^2)))#SE_A

  results[(smax-smin)/step+2,4]<-2*pnorm(abs(results[(smax-smin)/step+2,2]/results[(smax-smin)/step+2,3]),lower.tail=FALSE)#p-value for null=0

  results[(smax-smin)/step+2,6]<-sqrt(sum((y-results[(smax-smin)/step+2,2]-results[(smax-smin)/step+2,5]*x)^2)/((length(y)-dim(results)[1]+1)*sum((x-mean(x))^2)))#SE_B

  results[(smax-smin)/step+2,7]<-2*pnorm(abs((results[(smax-smin)/step+2,5]-1)/results[(smax-smin)/step+2,6]),lower.tail=FALSE)#p-value for null=1

  results[(smax-smin)/step+2,8]<-1-sum((y-results[(smax-smin)/step+2,2]-results[(smax-smin)/step+2,5]*x)^2)/sum(y^2)#R2

  results[(smax-smin)/step+2,9]<-((length(x)-2)/2)*((results[(smax-smin)/step+2,8]-(1-(sum((y-x)^2))/(sum((y-mean(y))^2))))/(1-results[(smax-smin)/step+2,8]))#F_test_R2_based

  results[(smax-smin)/step+2,10]<-pf(results[(smax-smin)/step+2,9],2,length(y)-2,lower.tail=FALSE)#F_test p_val

  return(results)

}

The Attempt at a Solution



I would like to obtain a MATLAB code.[/B]
 
Last edited by a moderator:
Physics news on Phys.org
I added code tags for readability -- please always use code tags to post programs here. Thanks.

So this is for schoolwork? Can you show us your Attempt at the Solution? We don't do your schoolwork for you here.
 
berkeman said:
I added code tags for readability -- please always use code tags to post programs here. Thanks.

So this is for schoolwork? Can you show us your Attempt at the Solution? We don't do your schoolwork for you here.
It's not a schoolwork.
I would like to perform an anlysis of two time series x and y, my aim is the perform a particular regression approach, on them.
 
Okay, let me move it to the Programming forum then...

ADD -- Have you started the MATLAB code? Can you post what you have done so far?
 
berkeman said:
Okay, let me move it to the Programming forum then...

ADD -- Have you started the MATLAB code? Can you post what you have done so far?
Not yet, i didn't understant what this codes do, i don't understand R programming language.
 
Tilfani said:
Not yet, i didn't understant what this codes do, i don't understand R programming language.
If you don't understand the R code, you will not be able to convert it to MATLAB code. If I were you I would start by looking at the R documentation.

The three lines below seem to be libraries that the code uses. If so, it would be useful to look at these libraries to see what the functions in them do, what parameters are required, and what kinds of values they return.
Code: 
require(tseries)
require(fracdiff)
require(matrixStats)

It would also be useful to understand what this notation means. I don't know R, so I can't offer any advice, but I would be looking to see what DCCA_beta_avg is, and what the notation <- means.
Code: 
DCCA_beta_avg<-function(y,x,smin,smax,step)
From the line below, it's possible that R uses <- for assignment operations, but I can't say for sure.
Code: 
XTX<-var(x)*(length(x)-1)
 
Step 1 is to find the person who wrote the R code without a single comment and poke him in the eye. Do you even know what the input and the output looks like? If so, you may be able to write MATLAB code without knowing every detail of the R code.
 
  • Like
Likes   Reactions: berkeman
Why do you want to convert it to matlab? It is possible to run MATLAB from inside R and use the results for computation. This would be probably easier.
 
  • Like
Likes   Reactions: berkeman
Mark44 said:
Code: 
require(tseries)
require(fracdiff)
require(matrixStats)
Yes, they can be found on the CRAN archive, including their descriptions.
 

Similar threads

MATLAB Problem with Runge-Kutta script file
  • · Replies 4 ·
Replies
4
Views
2K
MATLAB Transforming part of matlab code to Fortran90
  • · Replies 8 ·
Replies
8
Views
3K
MATLAB How do I write this chua oscillator to an audio file in matlab?
  • · Replies 5 ·
Replies
5
Views
3K
MATLAB MATLAB GPU Cross-Correlation Code Fails with OOM on 8 GB Supercomputer
  • · Replies 4 ·
Replies
4
Views
794
MATLAB How to plot flow data from an orifice in MATLAB without overlapping the plots?
  • · Replies 2 ·
Replies
2
Views
3K
MATLAB To what extent is the fit to experimental data good?
  • · Replies 12 ·
Replies
12
Views
4K
MATLAB Fixing a Matlab Code: Help Needed for Fourier Transform Calculation"
  • · Replies 29 ·
Replies
29
Views
5K
MATLAB What is the issue with my step function code and how can I fix it?
  • · Replies 3 ·
Replies
3
Views
2K
MATLAB Crank-Nicholson in 2D with MATLAB
  • · Replies 1 ·
Replies
1
Views
4K
MATLAB System of differential equations
  • · Replies 9 ·
Replies
9
Views
2K