Convert R Code to MATLAB for Detrended Cross Correlation Analysis

Click For Summary

Discussion Overview

The discussion revolves around converting a specific R code implementation for Detrended Cross-Correlation Analysis (DCCA) into MATLAB. Participants are focused on the technical aspects of the code, including the calculation of beta coefficients and various statistical measures derived from the analysis.

Discussion Character

  • Homework-related
  • Technical explanation
  • Mathematical reasoning

Main Points Raised

  • One participant presents the R code and requests assistance in translating it to MATLAB, highlighting the main functions involved in the DCCA process.
  • Some participants discuss the structure of the R functions, such as DCCA_beta_avg, DCCA_beta_sides, and DCCA_beta_SE, and their roles in calculating various statistical outputs.
  • Others propose potential MATLAB equivalents for specific R functions, considering differences in syntax and function availability between the two programming languages.
  • A later reply questions the handling of random number generation in the R code and how it might be implemented in MATLAB.
  • Some participants express uncertainty about the direct translation of certain statistical functions and suggest alternative approaches to achieve similar results in MATLAB.

Areas of Agreement / Disagreement

Participants generally agree on the need to convert the R code to MATLAB but express differing views on the best methods to achieve this. There is no consensus on specific translation strategies or the handling of certain statistical computations.

Contextual Notes

Limitations include potential differences in statistical function implementations between R and MATLAB, as well as assumptions made in the original R code that may not directly translate to MATLAB without modification.

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
4K
MATLAB MATLAB GPU Cross-Correlation Code Fails with OOM on 8 GB Supercomputer
  • · Replies 4 ·
Replies
4
Views
6K
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