Convert R Code to MATLAB for Detrended Cross Correlation Analysis

In summary, the conversation discusses converting R code into MATLAB code that performs a calculation of beta coefficients based on Detrended Cross Correlation Analysis. The code uses various functions and equations from the tseries, fracdiff, and matrixStats packages. The conversation also includes discussions on calculating standard error and p-values for the beta coefficients, as well as obtaining a weighted average for the coefficients.
  • #1
Tilfani
12
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
  • #2
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.
 
  • #3
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.
 
  • #4
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?
 
  • #5
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.
 
  • #6
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)
 
  • #7
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 berkeman
  • #8
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 berkeman
  • #9
Mark44 said:
Code: 
require(tseries)
require(fracdiff)
require(matrixStats)
Yes, they can be found on the CRAN archive, including their descriptions.
 
  • #10
  • Like
Likes FactChecker

1. What is Detrended Cross Correlation Analysis (DCCA)?

Detrended Cross Correlation Analysis (DCCA) is a statistical method used to measure the long-term correlations between two time series data sets. It is an extension of the traditional cross-correlation analysis, which assumes that both time series are stationary. DCCA takes into account the possibility of non-stationarity and trends in the data, making it a more robust method for measuring correlations.

2. Why would I want to convert R code to MATLAB for DCCA?

There are a few reasons why you may want to convert R code to MATLAB for DCCA. First, MATLAB has a more user-friendly interface and is better suited for scientific and engineering applications. Additionally, MATLAB has a larger library of statistical and data analysis functions, making it easier to perform complex analyses like DCCA. Some researchers may also have a preference for using MATLAB over R, so converting the code allows for a wider audience to use and replicate your analysis.

3. What are the main steps involved in converting R code to MATLAB for DCCA?

The main steps involved in converting R code to MATLAB for DCCA include understanding the code and its purpose, translating the code into MATLAB syntax, and debugging and testing the code to ensure it produces the same results as the original R code. It may also be helpful to consult MATLAB documentation or seek assistance from colleagues who are familiar with both languages.

4. Are there any limitations to converting R code to MATLAB for DCCA?

One limitation of converting R code to MATLAB for DCCA is that the languages have different syntax and functions, so not all code can be easily translated. Some functions in R may not have an equivalent in MATLAB, or the syntax may differ significantly. Additionally, the version of MATLAB being used may affect the compatibility and functionality of the converted code.

5. Is it possible to use both R and MATLAB for DCCA in the same analysis?

Yes, it is possible to use both R and MATLAB for DCCA in the same analysis. This can be useful if certain parts of the analysis are better suited for one language over the other. However, it is important to ensure that the results obtained from both languages are consistent and that any conversions between the two are accurate.

Similar threads

MATLAB How do I write this chua oscillator to an audio file in matlab?
    • MATLAB, Maple, Mathematica, LaTeX
Replies
5
Views
995
MATLAB Transforming part of matlab code to Fortran90
    • MATLAB, Maple, Mathematica, LaTeX
Replies
8
Views
2K
MATLAB How to plot flow data from an orifice in MATLAB without overlapping the plots?
    • MATLAB, Maple, Mathematica, LaTeX
Replies
2
Views
1K
MATLAB To what extent is the fit to experimental data good?
    • MATLAB, Maple, Mathematica, LaTeX
Replies
12
Views
3K
MATLAB Fixing a Matlab Code: Help Needed for Fourier Transform Calculation"
    • MATLAB, Maple, Mathematica, LaTeX
Replies
29
Views
4K
MATLAB Suggestions for a non-linear fit model in Matlab
    • MATLAB, Maple, Mathematica, LaTeX
Replies
1
Views
2K
MATLAB MATLAB Plug flow reactor optimisation Problem
    • MATLAB, Maple, Mathematica, LaTeX
Replies
18
Views
3K
MATLAB Writing mathematical equations in MATLAB
    • MATLAB, Maple, Mathematica, LaTeX
Replies
1
Views
3K
MATLAB What is the issue with my step function code and how can I fix it?
    • MATLAB, Maple, Mathematica, LaTeX
Replies
3
Views
2K
Mathematica Learning how to compute Christoffel symbols using Mathematica
    • MATLAB, Maple, Mathematica, LaTeX
Replies
6
Views
5K

Hot Threads

Recent Insights

Back
Top